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The work. described in this report was performed by the Control and. 
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ABSTRACT 

A detailed computer model for the rate kinetics of an atomic 
vapor laser excited by electrical discharge is proposed. The 
model equations are defined and the computer program structure is- 
discussed. 
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PROPOSED COMPUTER MODEL FOR ELECTRIC DISCHARGE 
ATOMIC VAPOR LASERS 


Introduction 

It is the goal of the present program of computer modeling to numerically 
simulate the gross features of atomic vapor lasers excited hy electrical discharge* 
Due to the inherent complexity of interactions in such a laser, the extent of 
attainment of this goal is uncertain* Nevertheless, much qualitative xnformation 
on some of the processes taking place can be expected, leading to a better under— 
standing of the physics of this laser type* 

Limitations on computation time and computer program storage make inevitable 
considerable simplification of real processes or characteristics of the laser* 

As a first step, the proposed model focuses on rate processes, i.e;, the change 
of electronic excited state populations through inelastic collisions and radia- 
tive interaction* Only temporal variations are included; the system is taken 
spatially homogeneous* In any particular atom or ion, levels of similar configu- 
ration (e*g* terms) are handled as a unit with populations assumed proportional 
to degeneracies* The set of level units of the atom or ion is split into three 
distinct groups: an upper subset in equilibrium with the free continuum states, 

a middle subset assumed quasistationary and a lower subset requiring numerical 
integration of rate equations to determine populations* The continuum plus the 
upper subset are treated together and will be called the ‘’extended continuum” [1] • 
The relaxation times between states in the extended continuum and also in the 
quasi stationary group are short relative to any other physical time of interest 
in the system* Populations in intermediate levels are small compared to the 
lower levels or the extended continuum. The occupation probability (based on 
equilibrium) is high for the low levels due to the Boltzmann factor and high for 
free or near free states due to large degeneracies. The short relaxation times 
and large possible population fluxes lead to equilibrium of extended continuum 
states among themselves* Fluxes in the quasistationary group are relatively 
small and the net rates into and out of a particular level are nearly the same; 
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populations change in a slow adiabatic manner so that they may be approximated by 
algebraically solving the rate equations assuming null rates [2], It is assumed 
in the model that the degree of ionization of the laser gas is always sufficient 
(greater than 10 to neglect the effects of atom-atom collisions on rates [1]; 
thus excited state equilibrium is characterized by the electron temperature 
without regard to heavy particle temperatures* This amount of ionization also 
assures a Maxwellian electron velocity distribution [2]. Neglect of atom-atom 
inelastic collisions in the rate equations does not mean that they are necessar- 
ily inconsequential in sublevel relaxation or relaxation between close levels in 
a unit* An equilibrium extended continuum also requires a sufficiently small 
charged particle mean free path. 


Rate Equations 


Let the population of level unit n, specie s with parent ion core charge 
ze be The extended continuum is denoted by n^; the corresponding reference 

^ 5 Z S Z’f'i 

number density is N ^ . Note that the actual density of the extended 

n^* 1 

continuum is 


(1 + ez 


3/2 


) N 


3 , Z+1 
1 


wnere 


e = N^(TTe^ « 1 [1], 


This holds for an extended continuum lower principal quantum number of about 5 
to 7 and is based on the Debye lowering of ionization potentials [3J. Here 
is the electron density and is the Debye length. A Saha factor Is defined as 


6 ®’^(T ) 
n e 


& 


s , z 
h 


2 g 


s , z+1 



(I®’^/kT ) 
n e 


where the g*s are degeneracies, h the Planck constant, the electron mass and 
is the ionization potential of level n. At equilibrium, = 6®*^ N N 

Q * n n u X 


- 2 - 



77-11 


with 0^!^ N E 1. For hydrogenic levels, the preceding degeneracy factor equals 
n '' e 

the square of the principal quantum number. If the ionization rate coefficient 

from n is (T ), the total excitation rate coefficient into the extended 

n» e 

continuum is 

xinx n«> X nm 
ra>n* 


Excitation coefficients K are discussed in Appendix D. The maximum principal 


nm 


1/2 


At equilibrium, N K, = N K , giving deexci ta- 
rn mn n nm 

E where A is the 
mn mn 


D- o 

The net radiative emission rate lii ^ n is A 


quantum number is (z5,^/2a^) 
tion coefficients. 

Einstein spontaneous emission coefficient and E is the escape factor (Appendix B) 
If a particular transition is coupled into the laser cavity , E is replaced by a 
factor dependent on the laser intensity (Appendix C) / ’In the following, E will 


denote a generalized emission factor. The emission rate depends on the transi- 
tion line profile unless the transition is optically thin; line profiles and 
broadening for the model are the topics of Appendix A. 


A norpia-lized population is defined hy ^ ^e^n*^ that at 

equilibrium (all li) . Also.R®’^ =(d/dt)F’^ and =(d/dt)5,n N^' 

are net rate functions . 


For n < n*, the basic rate equations for an atom or ion are 
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The quasistationary group is n* > n ^ n. The rate (d/dt)N^’^ =•• 0 for this group 

corresponds to setting R , R and d£nT /dt null in the above. Rate equations 

*n 6 s 

for all atoms and ions, plus those for cavity intensities (Appendix C) and the 
electron energy equation for dinT^/dt together define the model. An implicit 
assumption in the model is that for any given core, only one excited electron 
is in a discrete state, the rest being free. The radiative recombination term 
is given by H 1, 


h*n 


32 ct^Ryd ^4 
3 ^ 



where a is the fine structure constant, p^ the principal quantum number, the 
Gaunt factor, the energy gap to the extended continuum, and F(x) = e^Ej^(x) 

tl]. 


Charge and Specie Conservation 


The constraints produced by charge arid specie conservation are given next. 

It* is assumed that the highest stage of ionization of a specie atom is an ion 
ground state denoted by N^. Total specie number density is N and maximum parent 

g 

charge number is z , Given N , N is found 

S S X 


N = (1 + N® + N (s® + eS®) , 

s'l e\a b/’ 


< ^ L E . 


z-l n<n^ 


H x: <- - 1)'^' 8i’" 
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The terms proportional to e give the populations of the bound states parts of the 
extended continuum* The electron density is obtained using 


C - £A = N^( 1 + D - GB) , 


S 


s 

c 


Z E 1) • 

z n 




(z - 2) (z - 1) 


3/2 


gs,z -^,z 



B = 


S ' 


- z S® + z^'^^ 
d s b .s 




c = z K , 

^ s s * 

s 

D = E (z - sM 

\ s a c ^ 

c* ' 


In the above, G is calculated with N = C/(l + B) and then the final value of 

e ^ 

N is found* 
e 

Neglecting terms of order G, the conservation and rate equations may be ' 
combined to give 


\ = 


(Zs +1 


- z) 


(\ 


. s , z+1 s , z 

- "l “r 


)■ 


s,z 
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s,z ^ 

I n n nn 


z 
nn* 


is an ionization coefficient, and 


“r’“ - E ("eC' * «n’i 


n 


is a recombination coefficient. Similarly, ‘'rate heating" of the free electron 
eas D is given as 


°.-".E 4’“ - »e . 

s 


(n 

R / n n \ e nn* n* n*n/ n / 


n 


is a recombination heating coefficient and 


dS,z 


' E 4-' (4’^ ^ E <u^( 4’^ - 4-')) 


n 


k<n^^ 


is a coefficient of energy absorption for inelastic collisions. 
Reduction ,of Rate Equations 

The rate equations may be rewritten in the form (for given s,z) 

R = /.M,N, + F N + 
n nk k n 1 

k<n* 
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where the + superscript denotes z+1. Nulling and dilnT^/dt in matrix M gives 
matrix M\ The rate equations are of the stiff type, i.e,, characterized by 
multiple time scales. Integrators designed to handle this type of problem [4] 
require the Jacobian matrix J = 3R /3N to scale the numerical integration time 
step. For the quasistationary levels, = 0 for M replaced by giving 

k,m k 

where j,k ^ n and m < n. Substituting this expression for the quasistationary 
populations into the rate equations for the lower levels (n < n) gives the 
reduced set of equations 

R = / H N + G N,**" 

n ^ nm m n 1 

m<n 

where 

ri =- M - / ■ M '■ M* M,' - , 

nm nm n j j k km ^ 

j,k 

. G =. F - M . M' P, 

* n ’ n X ^ n J j k k 

The Jacobian of the reduced set P may be found with the aid' of the matrix , 
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Using the chain rule, 


P = J + 7 J . T. 

nm nm nj jm 

j 

and 


0 = J,’ + 

• km 


kj jm. 


(k ^ n) . 


(Primes have the same significance for J as M. ) 


Combining, 


P = J 
nm nm 


-E 


J . j! J' 

nj jk km 


where the matrix J is calculated for the quasistationary levels independent. In 
general, the reduction of the rate equations leads to two matrix inversions. 

The dependence of the rates on and leads to cross terms in the 
reduced Jacobian for different s,z. Due to the algebraic complexity of the 
model, the Jacobians are formed assuming fixed transition line profiles » The 
error in this assumption should not be large since the short time scales are 
associated with upper states which are usually optically thin. 


Energy Equation , and .the Penning Effect 


The electron energy equation is taken from 13-moment solutions to the 
Boltzmann* equation [1,5]. Included in the equation are the rate of change of 
enthalpy, elastic electron-heavy particle energy exchange, continuum radiation, 
collisional-radiative heating D^, and joule dissipation. Provision is made in 
the model to input either the current density or an imposed electric field or to 
bypass the energy equation by providing the temporal variation of the electron 
temperature. Heavy particle temperatures are assumed equal and constant. 
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Due to the possible importance of Penning excitation exchange in lasers 
having noble gases [6] , source code for one and two electron excitation and/or 
ionization of a receptor atom by a metastable donor is part of the model. This 
code is placed in blocks for easy removal. The probability of double ionization 
is apparently small [7]. 

Discussion 


The model for the laser rate kinetics has been designed to be of a very 
general nature. Tlie desire for flexibility and to test concepts led to the 
program form. Atom or ion level structure, oscillator strengths, atomic con- 
stants, etc. are input by distinct subprograms that may be changed at will to 
provide an arbitrary mix of elements and ionization stages. In order that this 
be done with maximum storage economy, this information is transferred to common 
blacks on an end-to-end basis. Further flexibility is obtained by the use of 
Univac Fortran V Parameter variables for dimension information, DO loop constants, 
index constants and the like. These variables are replaced by their assigned 
values at compilation. The specification statements for the Parameter variables, 
array dimensions, etc. are not part of subprogram source code but are placed in 
Fortran procedures (by the Procedure Definition Processor) for Inclusion at 
compilation. Thus only a single set of source code defining the procedures need 
be changed before compilation to modify the storage requirements of the assembled 
program. It is thus possible to minimise storage for a given set of atoms and 
ions with ease and without error, The numbers of level units in the divers 
groups can also be readily changed. 

The program is designed to treat lower level groups as stiff and Jacobian 
matrices are calculated. One purpose of the model is to determine if and/or 
when use of the quasistationary approximation for intermediate groups removes 
the stiffness so that simpler integrators may be used. Success of this tactic 
has been recently reported [8], but its use in cases of fast pumping is question- 
able [IJ. Use of the quasistationary group reduces storage requirements; it may 
or may not help insofar as reducing the program running time. It is necessary 
to determine the effect on running times and also constraints (c.g. values of 
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n and n*) of the three group model. Another item of consequence is the condi- 
tions. for an emission' factor E of unity . Calculations related to this factor 
are relatively extensive. After initial calculations, simplifications may ensue. 

A detailed listing of the model is provided in Appendix E. 

The program is in the process of being debugged. A sample run utilizing a 
single argon atom only with twelve level units has been apparently successful. 
Use of multiple species has led to difficulties with the matrix inversion aspect 
of the program; this anomaly and the configuration of output data remain as 
program problems. Numerical results will be given in a later report. 
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APPENDIX A 
LINE BROADENING 


It is desired to have a complete, yet compact, description of the excited 
states phase shifts (homogeneous broadening) and velocity changes (inhomogeneous 
broadening) » The most rigorous approach is that using the quantum density matrix 
[9,10]* The need for simplicity in the computer model leads to consideration of 
classical theory for the estimation of line profiles • 'In the model line center 
frequency shifts are neglected and the Voigt integral is used to enfold the 
Lorentz impact and Doppler profiles* The theory is presented in [11]? results of 
interest are summarized below* 

For an Inverse power law potential, V * emitter-perturber 

interaction the optical or phase shift cross section is where the Weisskopf 
radius 



The symbols are: e-^electron charge, a^ - Bohr radius, r - distance between 

interacting atoms, v - relative velocity of atoms, c - velocity of light, a - 

fine structure constant* Constant p is proportional to the critical phase shift 

m 

that defines the effective interaction range; = 0.318, = 0*411, pg = 0.533, 

Pl 2 = 0*773. This impact result is valid for angular frequencies in a region at 

the line center, Aco < v/b^, and for perturber densities smaller than 

These limits are of small consequence for the problem under discussion. The 

optical collision frequency (line semi-half width) is irN^^vb^^ where the velocity 

distribution is taksn as Maxwellian. For elastic collisions the phase shifts 

between the upper and lower states of a line subtract and C^. = - 

C ,, s . In the case of an ion emitter and electron perturbers, a correction 
ra( lower) 

for the hyperbolic collision path may be necessary [11]. If b^ is the Coulomb 
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collision radius, > b^, the cross section increases by a factor of approxi- 

g 

mately 2.3(b^/b^) , s = (2m - 4) /(2m - 3), giving an inverse square root 
dependence of the width on electron temperature. Interactions considered in 
the model are: 


1) resonance broadening (ground state atom perturber and like atom 

emitter), |c^j = (■nA'E,); (this interaction involves 
excitation exchange and is not elastic; the interaction of upper and 
lower states is added); 


2 ) 


electron-quadrapole broadening, |C 
*1^ L(L + 1)/(L(L + 1) - 3/4) for the emitter; 


“ A' 


where q 


2 

L 


3) , electron quadratic Stark broadening, 



sums are over all emitter bound states; 


4) 


5) 


Van der Waals broadening (ground state atom perturber and unlike 

atom emitter), C - a ; 

6 PVG//n ’ 


atom perturber and ion quadrapole emitter, 
(The polarization interaction C, 
and gives no broadening*) 




is state independent 


Broadening by, ions is small and is neglected. Symbols are: f , f , - 

ems abs 

emission and absorption oscillator strengths, Ryd - Rydberg unit of energy, AE - 
transition energy, g.,g - lower and upper degeneracies, average square of 

U N e/ 

the active electron radius of a state, L - total angular momentum quantum number, 
- perturber atom polarizability, Z - ion charge number. 
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The parameters in the resonance broadening constant refer to the 

transition between the emitter state and the (perturber) ground state. If the 

emitter state effective principal quantum number n > 3, the oscillator strength 

is small and the excitation exchange interaction is weak; further, the emitter 

radius is large. In this situation, repulsive electron exchange effects become 

1/4 

dominant [12]. For n (ac/v) , these effects can be estimated by using 

b. Effects at larger n [12] are not relevant to the model. The raaxi- 

W \e/ 

mum of the widths calculated using the two interactions is used. 


Broadening by electrons is estimated as the largest of the quadrapole and 
quadratic Stark interactions. Inelastic collisions between sublevels are fre- 
quently important in Stark broadening. Provision is made in the model to calcu- 
late the Stark width from curve fits to the calculations of Grlem [3], assuming 

k. IT 

additive upper and lower state interactions, Curve fits of the form CN^n T^ 
give good results for many levels (typically k = 5 and r = 0.4). 


The form of the Van der Waals constant G. is an approximation valid in the 

o 

limit of large energy level separations [3,11], As with the case of resonance 
broadening, electron exchange effects can be Important, especially for small a^. 
Only the simple form above is used; more precise analysis requires involved 
calculation and more knowledge of potentials than is commonly available. The 
polarizability a„ can be found from data [13] or estimated as 4aQ (Ryd/AE ) where 
the resonance level-ground energy gap pertains [3], Also, 


n^(5n^ + 1 - 35.(5. + 1))/2Z^ where Z refers to the parent ion. 


-I r^r-1 ^ 


Finally, the total homogeneous line width is the sum of the widths ofbtained 
from the formulas .above plus that due to the frequencies of quenching of the 
upper and lower states. The quenching mechanisms are that of spontaneous emis- 
sion (natural broadening) , induced emission or absorption (power broadening) and 
inelastic excitation and de-excitation between levels ("rate" broadening). Note 
that the quenching mechanisms give a Lorentz profile without constraint (this 
includes resonance broadening). 
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APPENDIX B 

RADIATIVE ESCAPE FACTORS 


An escape factor for a given transition gives the effect of induced radia- 
tive processes on the emission rate. The factors are calculated assuming a 
uniform medium and are averaged over the gas volume. For a non-inverted line of 
large optical depth, trapping gives an escape factor much less than unity. Con- 
versely, amplified spontaneous emission from an inverted line gives an escape 
factor greater than unity. 

The theory of line trapping is discussed in [14-19], Radiation loss from 
any particular line may be written as hv^ A^^ E where is the upper state 
density, hv^ the line center energy, the Einstein spontaneous emission coef- 

ficient and escape factor E =* <exp(-T^)>, being the optical depth (negative 
for ah inverted line), The average is taken over frequency and solid angle and 
for the model also over spatial position, A cylindrical geometry and a Doppler 
line profile gives E “ 0,90 /(t^ /in T^) where is the line center optical depth 
based on the radius, A Lorentz profile gives E == 0.63//F, Trapping decreases 
as drops to order unity and for small depth E =» 1 - In the model, the 

maximum of the Doppler and Lorentz factors (large limit, each calculated 
using the appropriate limiting line profile) is used to approximate trapping for 
the Voigt profile; then a smooth interpolation tO' the small depth limit is 
performed so that the entire range of optical parameters is covered in a con- 
venient yet reasonable manner. Note that at large depths, escape is determined 
by the line wings where the homogeneous broadening dominates since it is of 
distinct proportion. The approximation used gives neglible error at the larger 
optical depths where trapping is most pronounced, (Some error may occur due to 
non- impact optical interactions not properly accounted for in the model; the 
effect of this is considered slight since it can manifest itself only in cases 
of near complete trapping, i,e,, at small optical transition rates,) 

1 

An interpolation scheme between small and large gain limits is also used 
for inverted transitions, A discussion on the calculation of <exp(|T^|)> for 
It I > 1 follows. Let the cylinder center line optical depth value for the tran- 
sitioa line center be g, L the cylinder length (which defines g) , R the cylinder 
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radius, r the distance from an emission point to a cylinder boundary point, 
p(v) the line profile relative to the line center so that p(v^) = 1> the 

cylinder aspect ratio b =■ L/(2R), Then E - <exp(gpr/L)>. Averaging over fre- 
quency, E =<F(gr/L)> where 


F(s-) = 

/ TT 



(x) exp (s p(x)) dx 




k = /ir w p 5 w the Voigt half-width and p the line center absorption coeffi- 
v'^o V ° _ 2 -1 

dent, X = (V - V )/w . For a LorentJS profile p « (1 + x ) , 

F(s) = /ir k exp (s/2) I (s/2) k exp (s)/v^ for large s. For a Doppler pro- 

2 ^ 

file, p = exp(-x ) and it can be easily shown that F approaches the same limit 
for large s. Since F(0) = 1, the large gain limit is calculated using 
F(s) = 1 + k(e® - 1)/*^ In the Lorentz limit, this approximation gives a rela- 
tive error £20% for 0 < s £ 1, £5% for s > 2. In the Doppler limit, the error is 

as large as 34% at smaller s. The similar results for the profile limits is to 

be expected since it is the line center that controls the amplification process. 

Let the cylinder center line be parallel to the z-axis in rectangular coor- 
dinates, offset in the x direction by amount tsR, 0 £ Ci) £ 1, The equations for 
the cylinder wall are 

X = 2R q cos 0 

y = 2R q sin 0 

q E •- |u) cos 0 + (1 - sin^ 
z = r cos (}), 2Rq =* r sin (|) 

’ 7T 

Consider 0 £ z £ OL, 0 £ U £ 1; then for the cylinder side £ 4> £ where 
q ctn <t> = bU. The solid angle differential is sin ({)d<})d0. The emitter position 
is characterized by the variables u and u; the volume average is given by 
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integration over 2 tudcodu. For H(s) = (exp(s) - l)//s, the excess side loss' 


\s “ ~ ^^side aijproximately 


es 


("sin ♦ ) • 


The integrals over u) and 0 can be transformed to a single integral over q and 
the integrals over U and cj) can be transformed to a single integral by integrati 
by parts. Using u = gq/b, the result is 


-a ^Kb 
""es TTg 


j du A - (b„/g)2 J 




H(s) 


ds/ 

1 

A 

J 2 2 

s - u 


After removing the singularity at s = u by integration by parts, the remaining 
integration is done numerically. The results are given in Figure Bl, It is 
seen that is significant only for g > b, corresponding to large transverse 
gain. For the low gain limit ^es ^ ® 

last factor is 


- -2 
TTb 


f 


2 J 

q dq 


^ (t tan ^ I - i fn 


1 + 




After numerical evaluation for b > 1, the following curve fit gives this factor 
with an error less than 1%; 

<r/L>sidg = 0.66/b - 1.09/b^ + 2.30/b^ - 1.60/fa^ . 

The average solid angle that the side subtends is 4'IT where 


s 

TTb 


j dq A - q^ ( >/q^ + b^ - 


q) - 1 - 0.4244/b + 0.172/b' 


- 0.090/b^ + 0.025/b^. 
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If a is the ratio of the Lorentz half-width to Dopoler width, then 
k = (1 + 0.7304a + 0.5811a^)/(l + 1.2946a + 1.0299a^) 


and 

<^> i (0.707 + i. 25a)/(l + 2.5a) within an error of 2.3%. 
Similar calculations apply to the cylinder ends 5 



and 


<r/L> 


ends 


if 

✓ o 


qdq (cos ^ q - qN^l - q^) An 


1 + 


= (0.9428 + An b)/(8b^) + O.OlSA/b"^ 


Results of numerical calculation of the excess end loss are given in 

Figure B2. 
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APPENDIX C 
CAVITY EQUATIQNS 


Qaantum coherency effects are not considered in the modeling of the inter- 
action of strong cavity electromagnetic fields with atomic excited state popula- 
tions. Feedback from the optical cavity and saturation effects are estimated 
using rate type equations with steady state intensity profiles without velocity 
distribution or spatial hole burning. The present form of the cavity equations 
contain inherent inaccuracies that would be most pronounced for low pressure 
and/or very fast developing systems. This situation is dictated by the need for 
relative overall simplicity; future simplifications in other aspects of the 
model may lead to further development of the cavity equations. 

Denote the total intensity for a given* transition as I = chv N where N 

o ph ph 

is the photon density. The corresponding spectral quantities are I^ = chvN^. 

For each- direction of propagation, 

dN 3 


A is the Einstein emission coefficient, p the transition profile, r\ the solid 
angle and polarization factor, and the cavity decay time t^=L/(c £n(l/\/R^^) ) 
for length L and mirror reflectivities and Intensity is taken 

unpolarized and V replaced by except in the profile functions. (For plane 
polarized light, r| n/2.) Let the line center gain logarithm be 

Av the mode spacing, and the saturation parameter s E g/Zn (1/ Then 

summing over modes (frequencies) : 



(cos - 1) N / t + n A ^ 
ph c e uZ 


N 

u 
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where for mode intensities I 

m 


0) = 


2 ^ ' I p (V )/i p 
m in o 


and 


n 

, e 


23 nCVjjj)Av = /ndv 

m 


defines the total spontaneous emission noise into all active modes. Note that 


I = 


m 


and 03 < 1, The mode spacing Av = c£/(4L) where the cavity length and mirror 
radii define f. For the model, f = (x/(2 - x)) ' , 0 < x < 2, x is the ratio of 
length L to the effective mean mirror radius, (This equation defines the mean 
radius.) Only one mode exists if p^ Av > 1. The intensity profile- transition 
profile overlap Integral 05 is approximated by using steady state saturation 
intensity expressions, For homogeneous broadening a set of longitudinal modes 

gives [20]: 


03(s) 



As s ^ 1 (or p Av ■♦ “) , 0) ^ 1. This narrowing effect of saturation can be gener- 

O — 

alized by considering an intensity profile of the form p(v)/(l - sp(v)). Then in 

the limit s “^0, co <p> (a function of the ratio of Lorentz to Doppler widths. 

Appendix B) • Numerical integration is performed to find w in the limit 

p Av ^ 0, a curve fit to the results is modified for finite p Av to give 

(I3(s) = 1 - (1 - s)°‘^ (1 - <P>)/(1 + PqAv) . 
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This is the generalized saturation narrowing used in the model. 

The present cavity model can be expected to be a valid approximation at 
higher pressures where impact broadening is significant and at inversions not 
many times that at threshold. Different transverse modes have different diffrac- 
tion losses; the number of modes taken must be equal to those of relatively small 
loss. This number would define the effective emission solid angle It must 

depend on the number of photon passes in the cavity j cavity geometry, wavelength, 

and transition width. At present the parameter n is input into the program. 

2 ^ 

The single mode minimum value is (X/(2 ttR)) ; an upper bound is the average end 

solid angle (Appendix B) » Spontaneous emission noise is most important at the 

start of the intensity buildup and the larger values of ri are most pertinent. 

e 

The model presents the maximal effect of saturation narrowing within the 
rate approximation} calculations with oj(s = 0) may be used to determine the 
minimal effect. 
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APPENDIX D 

CLASSICAL EXCITATION RATES 


Expressions are needed for excitation rates between many pairs of levels. 

The expressions must be general in nature and relatively simple, yet reasonably 
accurate. The use of classical theory 'produces expressions that satisfy these 
constraints [1], Provision is made in the model to use empirical excitation rates 
from ground states. 

The present model utilizes the cross sections of Gryzinski [21]* Excitation 
from a level "A" to all levels "u" and above (Including the continuum states) is 
given by the coefficient [22] J 

where g^ is the number of equivalent electrons in the. lower state, is the 

energy gap, y 5 AE^^/kT^, A S 1 ^ 1 ^ 1 for lower state ionization energy- I^. 

The rate coefficient to level "u** alone is obtained by subtracting a similar 
expression with u -»■ u + Ij ~ \u-fl' Siven by K with 

A = 1. For A < 10 and 0.01 < y < 10, a good analytical approximation to numerical 
evaluation of rate coefficients using a Maxwellian velocity distribution is [22]: 


R = (3.84 X lO"^ cm^ eV^^^ sec"^) y*^ e"^ Cy^ + 7y/4 + 1/9)) , 

t = (A + 30)/(10A + 25). 

Asymptotic expansions for small and large y yield respectively 

R (y small) = 4.39 x i0“^ n/T (l - A In (Ay)/3) , 

R (y large) = 2.93 x lO"^ e“^ / (vTTy^), S = (3A -b 1)/(2A + 2). 

These functions are suitably matched in the intermediate regions. The above 
excitation function may be applied to allowed transitions. Gryzinski's exchange 
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cross sections give coefficients for forbidden transitions; a curve fit is made 
to numerical calculations [23]# 


Another classical calculation is that of Mansbach and Keck [24], where 
Monte Carlo trajectories are used in a three body (free electron, active electron, 
parent ion) system# For principal quantum numbers n^ and n^, excitation rate 
coefficient 


%u = 3.73 X lo“^ sec' 


■‘(r) 


0.17 n 


4.66 


H 


,-y 


n 


and for ionization 


^£oo = 1.87 X 10 ^ cm^ sec 


-1 I e_\ 

iRydl 


4.66 -Ay 
e 


(Calculations are made for hydrogen-like atoms.) These rates are smaller than 
that of Gryzinski for small energy gaps, AE ^ kT^ [24]. Since coefficients are 
large under this condition and the levels involved tend to be quasi-stationary 
[1], the difference is probably of small significance# Within the classical 
approximation, the Mansbach-Keck rates may be regarded as theoretically more 
rigorous. However, collective interactions between free electrons and a highly 
excited active electron make the accuracy of any lone perturber collision theory 
doubtful for transitions between levels of large quantum number# The model may 
be easily modified to employ Mansbach-Keck or quantum Born excitation rates. 


The form of the excitation coefficient as given by the R-function can be 
expected to be approximately correct for ions as well as atoms [25]. 
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APPENDIX E 
PROGRAM LISTING ' 

The computer program is written in terms of the following units or reference 
values for physical quantities. 

pressure 
density 

rate coefficient 
line width, relaxation time 
temperature 
atomic mass 

length (except wavelength) 
wavelength 

elastic cross section 
current density 
electric conductivity 
electric field 
magnetic field 
intensity 

The following symbol conventions are used generally throughout the program. 
The few exceptions are noted later. 

S 2 

AN vector of ground state populations of atom or ion type s,z 

»■ g 2J 

ANB vector of normal excited state populations N ’ 

n 

ANC vector of specie densities N 

s 

ANE electron density N 

® s 

ANP vector of densities of maximum ionization Nj^ 

CAR laser cavity aspect ratio b 

CL cavity cylinder length 

CMR cavity mean mirror reflectivity vTlj^R^ 

NEAR vector of level' indices n 

NCM maximum number of species (distinct nuclear cores) 

NSP number of atom and ion. types 
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NSTAR 

T 

TE 


Z or ZC 
ZS 


vector of level indices n* 
heavy particle temperature 
electron temperature 
charge number z of parent ion 
maximum charge number of specie 


In the Fortran procedure source deck SPECS, values of Parameter variables 
are to satisfy the constraints: 


If 


NLj^ ’= maximum number of level units of atom or ion type k as given by the 
data input subprograms, 

KM = maximum number of equations to be integrated, 


Then 


NLM 2 max (NL^) , 


ND > 



NL, 

k 


» 


M > 1 + NSP + 



+ 


3) 


2 


NT > max ((NLM + 3)^, NLM^ + NLM + 2o) , 

• KSV > 1 + 3 NLM + 5 NLM^/4 . 

Cn the format specifications in procedure SPEC3, the second record of format 62 
md the third record of format 63 refer to the vector ANC and must have data 
ength equal to NCM; the third and fourth record of 62 and the fourth record of 
>3 refer to vectors NEAR, NSTAR with dsta lengths determined by NSP. The program 
-S written for a maximum of four cavity equations. This is reflected in the 
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last dimension of ICAV, the dimension of WL, the last record of format 62 and 
corresponding last record of 63 referring to ICAV, the last dimensions of GP and 
LU in SPEC4 and the dimension of ANPH in SPEC7. Symbols in the procedures are 
defined in this appendix as the need arises. Common blocks LEVELS, EESND and 
PENCOM are used for input data and calculated primary variables, TEMPS, INCR, 

STORE and MSTOR are used for scratch and data transfer between subroutines. 

Along with the proper specification statements, .a particular problem requires 
that appropriate data and function subroutines be included in the program assem- 
bly. These subroutines are described in succeeding pages. Their specific names 
are required by the Collector for substitution in dummy subroutine calls. The 
proper Collector directives are essential to a problem. 
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-P0P»IFL SPECi> 

SPECIPROC 

PARAMETER RSP=3 »M=^1591»ND= 60*N».M=20 * NT=529 .NPL=NLM+ 1 ♦NSG=NLM»NLM» 
C MB=4*ND 

DIMENSION DS(M) *IND3(MJ , INX( NSP » 8 ) »BP ( MB ) »TMP(MT) 

COMMON/LE VELS/ 1, NX » DS , BP /TEHPS/TMP 
EQUIVALENCE (DS»INDS) 

END 

SPEC2 PROC 

PARAMETER NCM=3 

DIMENSION IC(NCM) .ANCCNCM) »ANB( ND ) »NBAR ( NSP ) »NSTAR(NSP) »ICAV{3*4)i 
C CFST(4) *3T(NDJ *rtL(A) 

COMMON/EESND/ICtNC»NBAR>NSTAR,ICAV»NCAV»CL»CAR»CMR»ETA,FCAV»CGM. 

• C PS»CFST*TE*T»ANE*ANC»ANB,BT,WL 
END 

SPEC3 PROC 

62 FORMAT ( 1P8E 10.3/ 1P3E10. 3/31 3/31 3/121 3 ) 

63 FORMAT UHl,2XlAhCAVlTY LENGTH=1PE10. 3 »2X16HLENGTH/DI AMETER= 

C lPElU.3»2Xl3riMIRROR REFL.=1PE10. 3 t2XllHM0DE RATIO=1PE10*3 » 

C 2X11HRES.. RATIO=1FE10.3/ 

C 3X3HTE=1PE10.3»2X2HT=1PE10.3»2X3HNE=1PE10.3/ 

C 3X3HNC=3{2X1PE10.3)/ 

C 3X5HNBAR=3(2XI3) ,2X6riNoTAR=3(2XI3 )/ ‘ 

C 3X16HCAVITY INDICES- 4(2x313) I 

END 

SPEC4 PROC 

PARAMETER NHSQ= ( NLM+NSQ ) /2 »NP4= NHSQ+2 1 

DIMENSION AN(NSP) *n(NHSQ) »WD( NHSQ ) »PH I ( NHS3 ) ,R ( ND ) , 

C GP(3.4) »LU12»4) 

CO^^MON/INCR/R 

EQUIVALENCE (WD ( 1) ♦TMP(2l) ) » (WD(1 ) »PHI { 1) ) , (W( 1 ) »TMP(NP4) ) , 

C (GPa»l)»TMP(l)),{LU(l>l)»TMP(13) ) 

IFN( 1»L)={ ( I-1)*(2*L-I) )/2 

END 

SPECS PROC 

logical lpn 

COMMON /PENCOM/ LPN » RPEN » JRC » JRL , J DC » J DL 

END 

SPEC6 PROC 

PARAMETER KSV= 561 
DIMENSION SV(KSV) 

common/store/sv 

END 

SPEC? PROC 

PARAMETER KM= 50 » NRSTM= 2*-NLM*NLM+KM*KM + 3 » ( NLM+KM ) 

DIMENSION ANR(KM) »ANPH(4) >STRMT X ( MRSTM) 

C0MM0N/MST0R/STRH7X 

END 
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Subroutine DATAIN is used to input and collect the data defining the particu- 
lar problem at hand. The. dummy calls LDUMYn are made equivalent to the names of 
the atom and ion data subroutines used in the problem by the Collector processor 
which generates the executable absolute program. It is assumed NSP ^ 6. The 
argument list in the data subroutines is as follows: 


■PQ 


AQ 

AQT 

DG 

ID 


W 

FR 

NL 

GNLI 

GFF 


vector effective principal quantum numbers based on term 
values 

vector angular momentum quantum numbers (active electron) 
vector total angular momentum quantum numbers 
vector degeneracies 

integer vector with the properties; ID(1) is the number of 
ground level equivalent electrons; otherwise, nominal 
ID(j) = 1, ID(j) < 0 if the Stark broadening of level j is 
to be calculated by the curve fit to Griem's tables (Appen- 
dix A), |lD(j)j = 2 if empirical excitation rates from the 
the ground level to j are to be used 
vector Gaunt factor for radiative recombination into 
levels [l] 

array of transition oscillator strengths, these are input 
^ 0 for nonallowed transitions (array is two dimensional) 
dimension of arrays, i.e., maximum number of levels 
degeneracy of ground level of next ionization stage 
overall Gaunt factor for recombination radiation. 


5 

6 


j<n* 


-3 

w.G.p . 

2 J j 


for pj = PQ(3), 


(I) is the relative weight of level j in its shell-, and for 

j 

hydrogenic levels the shell Gaunt factor 


G. = 1 - o.nspT^*^^ - o.osop"^^^ 
2 2 2 
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IS 

AM 

ALPHA 
CS, RS, KS 


integer parent ion charge number 
atomic mass 

polarizibility (applies to atoms only, null for ions) 
Stark curve fit parameters (Appendix A) 


The calls place data into temporary storage TMP, then a call to TRANS shifts the 
data into the array DS (INDS). DS(1) = NSP, DS(k + 1) = NLj^ (1 1 k £ NSP). 

The array INX(k, i) contains storage locations of the first elements of arrays 
for type k according to i=l: array is (IS, AM, GNLI, GFF, ALPHA, CS, RS, KS), 

i=2: PQ, 1=3: AQ, i=4: AQT, i=5: DG, i=6: ID, i=7: W and i=8: FR, The 

oscillator strength for a transition ra ->■ n 

|FR(m,n) | = |dS (m - 1 + NLj^ » (n - 1) + INX(k, '8)) | . 

For level units made from levels of different momentum quantum numbers S., the 

average of the product 8.(5. + 1) is used to determine AQ(j) For each type k, a 

set of four vectors of broadening parameters is placed sequentially into the 

array BP. The first vector is the resonance broadening constant C , the second 
2 2 

vector <r >/a , the third vector Stark constant C, and the fourth vector the sum 
e o’ 4 

of A-values to lower levels, i,e,,‘ the natural line width (Appendix A), Func- 
tions FI and F2 are used to complete the sums over levels in for j ^ n*. 

These functions are based on hydrogenic oscillator strengths. The vector IC 
gives the type index for atoms of the species and NC is the number of species. 

Array ICAV inputs the transitions that are to be considered in the cavity 
equations. The' first index of ICAV indicates the following; I for type index 
(or null) , 2 for lower level index and 3 for upper level index. NCAV is the 
number of cavity equations. Symbol ETA is the effective noise emission factor 
for the cavity and CGM is the ratio of the length to effective mean mirror 
radius (Appendix C) . The subroutine changes CGM to the mode spacing and calcu- 
lates the cavity decay rate as FCAV, the number of photon passes as PS and the 
wavelengths as WL. 

The dummy call to SDUMY is to be replaced by the Collector by EQUIL or QSTAT 
which initialize the level populations according to equilibrium or quasistationary 
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Q^j_ciil&tions iTBspscfcivBly » Thss© subjrovitinBS bitb dlsciissBd tli©y plscB 

the Saha factors in the vector BT. 

The Penning effect inputs into PENCOM are logical constant LPN — true if 
the effect is to be considered, RPEN - Penning rate constant (negative for two 
electron excitation), JRC - receptor specie index, JRL - receptor level index, 
JDC — donor specie index and JDL - donor level index. The donor specie must 
follow the receptor specie at data input. 

Subroutine DATAIN and corollary subroutines may be placed in a separate 
segment in the absolute program since they are needed only at the program 
start. Thus they do not represent any added storage in the program. 
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-TFOR»SI SUbl 

SUBROUTINE DATAiN 
INCLUDE SPEC1 *lIST 
INCLUDE SPEC2»LIST 
INCLUDE SP£C3»LIST 
INCLUDE SPEC5»LIST 
INCLUDE SPEC7*LIST 

PARAMETER NiSi = NPL+9,NlS2 = NiSl+NLM»N.lS3=NlS2+NLM,NlSA = NlS3 fNLM, 

C N 1 S S = N 1.S A+N LM » M S.6 = N 1 S 5+NLM 

DIMENSION PU(NPL) *AQ(NLM) *AQT(NuM) »D6(NLM) » ID(NLM) »W{.NlM) »FR(NSQ) 
DIMENSION R(ND) 

COMMON/ INCR/R 

EQUIVALENCE (IS* TMP ( 1 ) ) ♦ ( AM,TMP ( 2 ) ) ♦ (GNLI ♦ TMP ( 3 ) ) * ( GFF » TMP { A )) > 

C (ALPHA»TMP( 5) ) » (CS»TMP(6) ) » ( RS»TMP (7 ) ) > ( KS»TMP ( 8 ) ) ♦ (PQ( 1 ) » 

C TMP(y n *(AQ (1) »lMP(NlSi) ) ,( AQT(1 ) ,TMF (N152) ) » ( DG ( 1 ) » 

C TMP(N1S3) ) » ( ID(1) »TMP(N1SA) ) ♦ {W( i } ,TMP(N1£E ) ) » 

C (FR(1) ,TMP(N1S6) ) » ( JS*NSTRT ) , ( ANR ( 1 ) »£TRMTX ( 1 1 ) » 

C (ANPHiDjSTRMTXU) ) 

110 FORMAT(L7*iPE10.3,Al3) 

111 FORMAT ( 1HU,2X23HPENNING KATE CONSTANT = IPEIO. 3 > 5 XIOHRCPT CORE 13* 

C 5X11HRCPT level I3*5X9HDiMR CORE I3»3X10HDNR LEVEL 13) 

F1(X»Y)=1.0/(1.0-X*X/(Y#( Y-1.0) ) )«*A-1,0 
F2(X*Y)=0.A9*X**3*(i,U-2.8*EXF( -1 .A* AMAXl ( Y-X » 0. 75 ) ) )« 

C Fi{X»X+AMAXl( Y-X*o*75) ) 

INDSd )=NSP 
NSTRT=NoP+2 
L=1 

C DUMMY SUdR CALLS ARE TO BE CONVERTED BY COLLECTOR 

CALL LDUMYl(P(d»AQ»A(wT »DG» iD*N»FR»NL*GNLl ♦GFF»IS,AM»ALPHA*CS*RS,KS) 
IF(L*GT.NSP.OR.NSTRT.GT«M:-a-6*Ni.-NL*NL) GO TO AO 
CALL TRANS 


5 


IF(NSP.EQ.l) GO TO 5 

CALL LDUMY2{PQ*AQ»AQT>DG»ID»W»FR»NL»GNLI »GFF*IS»AM,ALPHA>C5»RS*KS) 
IF (L.G r.NSP.OR.NSTRT,GT.M-8-6*NL-NL«NL) GO TO AO 
CALL TRANS 

IF(NSP.EQ*2) GO 10 5 

CALL LDoMY3(PQ*AQ»AGT»DG* ID»W»FR»NL»6NLI »GFF * I S , AM, ALPHA ♦ CS *RS *KS > 
IF(L.GT.NSP.0R*NSTRT*G1 .M-8-6*Nt_-NL«NL) GO TO AO 
CALL TKANo 

IF(NSP«EQ*3) GO TO 5 

CALL lDUMYA(PO,Au»AQT»DG, ID,W,FR,NL*GNLI * GFF , I S , AM , ALPHA , C3 * RS , KS ) 
IF{L*GT.NSP.CR.NSTRT.GT*M-8-6»NL-NL»NL) GO TO AO 
CALL TRANS 

IF(NSP*EQ,A) GO TO 5 

CALL LDUMY5 (PQ,AQ,AQT,DG* ID*W,FR*NL*GKLI »GFF ♦ I S , AM, ALPHA ,CS ,RS,KS ) 
IF(L.GT.NSP.0R.NSTRT,GT.M-8-6*Ni-“NL«NL ) GO TO AO 
CALL TRANS 

IF(NSP*EQ.5) GO TO 5 - 

CALL LDUMY6(PQ,AQ,A3T,DG, I0*W*FR,NL*GNLI *GFF , IS,AM,ALPHA,CS*RS,KS) 
rF(L*GT*N3P.0R.NSTRT*GT.M-8-6«Nj--ML*NL) GO TO AO 
CALL TRANS 
JS=0 

DO 30 K. = 1*NSP 
L=INDS(X+1) 

L1=INX(K*3)-1 
.L2 = INX(K,2)-1 
L3=INX(K*5) 


LA=INX(K»8)-L 

Z=lNDS(L2-7)**2 


EEPROBUCiBILEPy Wdl 
ORIGINAL^:pA^B‘jjS!EOQB 



77-11 


DO 20 1=1, L 

, A= Z/DS{L2+1)**2- Z/DSC I+L2)**2 
IFU.EQ*!) A=1.0E9 
B=SQRT(DS(l3>/DS( I+L3-1) ) 

BP{ I+JSI = i,91*3*ABSlDS[ I*L+U) ) /A 
A=DS( I+L2 )**2 

B=1.0-3.0*DS( I + Ll)*{1.0+DSn+Ll ) ) 

BP ( I +JS+L ) =A* ( 5 • Q*A+B ) / ( 2 . 0*Z ) 

B=0,U 

C=0.O 

DO 10 J=1,L 

A= Z/DS( J+L2)**2- Z/DS( I+L2)^^*2 
IFU.EO.J) A=1,0E5 
A=SI6N(A*A,A> 

IF(l.GT.J) C=C+A*AB3(DS( I+L^^J+L4-1 ) ) 

10 B=B+2.o*AbS(DS( 1+L*J+L4-1) )/A 
B=B“F2(DS( 1+L2) tD3( LI) ) 

BP( I+JS+2*L)=B 
20 BP( I+JS+3*^L)=6.39*C 
JS=JS+4*L 
30 CONTINUE 
GO TO 50 
40 WRITE(6»42) 

42 FORMAT ( lHo,3X20riliMCORRECT INPUT DATA) 

STOP 
50 JS=o 

DO 60 K=1,NSP 

L=1NX(K,D 

L=1NDS(L) 

IF(L.Nti.l) GO TO 60 
JS=JS+1 

IF(-JS.GT,NCM) GO TO 40 
IC( JS)=K. 

60 CONTINUE 
NC=JS 

RE AD { 5 , 62 ) CL » CAR , CKR , E T A , CGM , T E , 1 , A NE , ANC , NBAK , NSTAR » I CAV 
DO 66 K=1»ND 
66 R{IO=0.u 
NCAV=0 
DO 70 K=l,4 

7u IF( ICAV( l,K).GT.O) NCAV=NCA\/+1 
. , DO 80 K=1,NSP 

IF(NSTAR(K)*GT. (INDS(K+1)+1) ) NSTAR ( K )= INDS ( K+1 ) +1 
8u IF(NbAR(K) .GT.NSTAR(K) ) NBAR( K.) =NSTAR ( K ) 

DUMMY SUdR calls ARE TO BE CONVERTED BY COLLECTOR 
CALL SDJMY 

WR I TE ( 6 , 63 ) Cl » CAR , CMR, ET A ,CGM, TE , T , ANE , ANC , N3AR , NSTAR , I CAV 
A=TVEF(1.0,CAR,1.0) 

IF(AbS(CGM-1.0) *GE.1.0.0R.CMR.GT.1.0.CR.CMR.LE.O,0) GO TO 40 
CFST{1)=CSA(CAR»0) 

CFST(2)=CSA(CAR,1) 

CFST(3)=CSA(CAR»2) 

CFST(4)=CSA(CAR»3) 

READ ( 5 » 1 1 0 ) LPN » RPEN , JRC » JRL , JOC ♦ JDL 
IF(.NOT.LPN) GO TO 96 
IF(ABS(RPEN).GT.0*0) GO TO 90 
LPN=*FALSE. 

• GO TO 96 
90 K=IC(JRC) 

- 35 - 



92 

94 


96 


97 

98 

lOU 


4 


11 


21 
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IF(RPEN*LT,0.0> go To 92 

IF( JRL.6T.NSTAR (lU ) JRk»NSTAR(Kl 

GO TO 94 

IF(K.EQ.NSP) GO TO 94 
IF(JRL*GT*MSTAR(K+1) ) JRL=NSTAR(K+1) 
WRITE(6»111) RPEN>JRCtJRl.»JDC»JDL 
IF( JRC.GE..JDC) GO TO 40 
IF( JDC.GT.NC) GO TO 40 
KrlC( JDC) 

IF-( JOL.GE.NSTAR(K) ) GO TO 40 
FCAV=299.7925»ALOG(1.0/CMR)/CL 
CGM=SQRT(CGM/(2.0-CGM) ) 
CGM=74.948»CGM/CL 
PS=1.0/AL06{1.0/CMR) 

IF{NCAV.LE.0) go to 98 
DO 97 J=1,NCAV 
K=ICAV( i*J) 

L=INX(K»1) - 
C=iNDS(U**2 
L=INX(K*2)-1 
I=ICAV(2»J) 

K=1CAV(3,J) 

A=C/DS( I+L)**2“C/DS(K+L)**2 

WL( J)=9.11267E-2/A 

L=0 

DO 100 X=1,NSP 

L=L+NBAR(K)-1 

L=L+NCAV+1 

IF(L.GT.KM) GO TO 4C 
RETURN 

SUBROUTINE TRANS 
INXU.*1)=NSTRT 
INX(L»2)=NSTRT+8 
INX(L,3)=NSTRT+NL+9 
DO 4 K=4,8 

INX(L»K)=INX{L»K-1)+NL 

L=L+1 


INDS(L)=NL 

INDS(NSTRT J =IS 

DO 11 K=2»7 

DS{K+NSTRT~i)=TMP(K) 

INDS{NSTRT+7)=KS 

OS ( NSTR T+NL+8 ) =PQ ( NL+ 1 ) 

DO 21 K=1»NL 


.pools. 


DS(K+NSTRT+7)=PQ(K) 

OS ( K+NSTRT+NL+8 ) -AQ ( K ) 
DS(K+NSTRT+2#NL+6 )=AQT ( <) 

DS ( K+‘NSTRT+3*NL+8 ) =DG ( A ) 

INDS ( K+NSTRT+4*NL+8 ) = I D ( K ) 
DS{K+NSTRT+3*NL+8)=w{K) 

DO 21 J=1>NL 

DS( J+NL»K+NSTRT+5*NL+8) =FR(J+NL*K-NL ) 

NSTRT=NSTRT+(NL+3)**2 

RETURN 

END 
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Function ERFGS collects empirical ground state excitation rates for atom 
or ion types k. These rates are placed in function routines and depend on the 
electron temperature and level index n. The functions are to return a negative 
value (e.g. , -1.0) for levels whose rates are not calculated. The Collector 
replaces the dummy names EDUMYi by the functions which are associated with the 
replacements for LDUMYi. 

Function ERFGRY is used to calculate excitation coefficients based on the 
classical Gryzinski cross sections (Appendix D) . For nonallowed transitions, 
function EREG based on Gryzinski exchange cross sections is used. 
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-TFOR*SI SUB2 ' i 

>- FUNCTION ERFGS(TE.N.K) 

L J . 6 0 TO g • 2.» 3i4=i5 » 6 ) » K 

’ C DUMMY FUNCTION REFS ARE TO BE CONVERTED BY COLLECTOR 
* 1 ERFGS=EDUMY1(N»TE) 

R^TJJRN _ 

2. ERFGS»EDUMY2(N*TE) 

RETURN 

. 3_ ERFGS»EDU MY3(N»TE) 

RETURN 

4 ERFGS*EDUMY4(N»TE) 

RETURN . . 

5. ERFGS=E0UMY3(N*TE) 

RETURN 

.. 6. _ERF0S=EDUMY6.(N»TE) 

RETURN 

END 
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-TFORtSi SUB3 

FUNCTION ERFGRY(Y»A) 

REAL Y»A 
IND=0 

IF(Y.LE.1.0E-3) GO TO 100 
IF(Y.GE.3C,0) GO TO 3C0 
IF(Y.GE.O.l.AND.Y.LE.lO.O) GO TO 200 
IF{Y*6T«lU.u) GO TO 5o 
X=ALOG10( Y) 

IF(X«LE.-0*iJ-0.25*A*AND.X,LE.“2.0) GO TO 100 
IF(X.GE* U.5-0.25*A,AND,X.GE.-2.0) 60 TO 20C 
IND=1 

IF(A.LE.6*0) GO TO 10 
IF(A«6E*10.0> GO TO 20 
WS=X+u.3+u.25*A 
GO TO lOO 
lO WS=X+2.0 
GO TO lOO 
20 WS=X+3.0 

GO TO Iv^O 
50 IND=2 

WS=o.05*Y-0.5 
GO TO> 200 

100 F=SQRT{ Y)*( 1.0-A*ALOG( A*Y) /3.0) 

IF(F.LT.O.O)- F=0.0 
IFdND.EQ.O) 60 TO 5‘00 
FT=F 

200 T=(A+3o,0>/(10.0*A+25.0) 

F= A**0 • 25* ( 1 . lA3*Y*Y+2 , 0*Y+0 .127) 

F=Y**T/F 

IF(IND*EQ.O) • GO TO 500 
IF(IND.E0.2) GO TO 250 
210 F=WS*(F-FT)+FT 
GO TO 500 
250 FT=F 

300 T=1.5-1.0/( A+1.0) 

F=1.5*SQRT(A)*Y**T 

F=1.0/F 

IF(IND.EQ.O) GO TO 500- 
GO TO 210 

500 ERFGRY= 0.875*F*EXP (- Y ) 

RETURN 

END 
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-TFOR»SI SUB4 

FUNCTION'EREG( Y*A»b) 

X=A/10u.O 

IF(Y.GT.29*0) GO fO lU 
C=0. 8665+ ( 4*634-2. B34*Y )*Y^^EXP( -3.0* Y )+l*288*Y**3*EXP( -Y) 
G=5.97*EXP(-Y)-10.65*EXP(-2.0*Y ) +7.8 5*EXP ( -3 . 0*Y ) ! 

GO TO Zv 
10 C=0,8665 

G=0.0 

20 X=0.524/X-2. 195+2. 493*X-0.397*XJ^X 

W=2.0*C+2.0E-2*G*X 
X=8/10.0 
FX=X-0,1 

FX=6.5*FX+lU.0*FX*FX 
I.F(FX.GT.88.0) FX=88.0 

C =(252. 26+3. 713/X+0.0307/X**2)*(1.0-EXP(-FX) ) 

IF(X.GT.3.1> GO TO 30 

G =4.U+0.0541*X+0.247*X*X»EXP(-4.0*X*X}+184.0*X**4*EXP(~9.C*X*X) 
GO TO 40 

30 G =4.0+U,U541*X 
40 FX=2. 013+0. 2098*X-0.0802*X*X 

G = G*(Y-1.0) 

IF( G.GT.88.0) G=88.0 

EREG=2.0E-3* C*Y**FX*EXP(-G )/A**W 
RETURN 
END 
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Subroutine NWTH calculates the line width per unit perturber number density 
for atom perturbers* The perturber type xndex xs KP and the radiator index is 
KR, NL denotes the number of radiator levels calculated. If the perturber is an 
ion, the subroutine sets logical constant TEST to true and returns. The width 
of the transition between lower level 1 and upper level j is placed in the vector 
W residing in TEMPS common. The storage is arranged in terms of a half-matrix 
(null below the diagonal). Subroutine EWTH performs a similar calculation for 
electron perturbers. The basic broadening theory is outlined in. Appendix A. The 
results of these subroutines along with natural, power and rate broadening are 
combined in LINWID. Power and rate quenching frequencies are obtained from INCR 
common. The total homogeneous broadening W and Doppler broadening WD reside in 
TEMPS. 

Subroutine CRATE calculates the collisional rate coefficient matrix 
for type KR. The matrix resides in TEMPS as the symbol RM. It is assumed that 
the equivalent number of ground state electrons (ANG) is contained in empirical 
excitation rates. Ionization rates are the matrix elements with j = NPL. 
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-tforjsi subs 

SUBROUT I NE NWTH ( T * KR , KP * NL » TEST , 

INCLUDE SPEC1,LIST 
PARAMETER NHSQ= ( NLM+NSQ ) /2 
dimension W(NHSQ) 

EQUIVALENCE (IS »TMP(1) ) » ( L»TMP( 2 ) ) . ( AMR» TMP ( 3 ) ) ♦ ( AMP »TMP ( A ) )» 

C (ALPH»TMP(5 ) ) . IRM»TMP(6) ) *(K»TMP (7) ) » (C6»TMP(8) ) »( C5P » TMP ( 9 ) ) » 
C (C3*TMP(10) ) » (C3'P»1MP(11) ) » (CU*TMP( 12.) ) »-{-CUP-,TMP {13 ) ) ♦ 

C (W( 1) »TMP(21) ) 

LOGICAL TEST 

F(X)sX*(X+1.0)/ (7.5*X*(X+1.0)-5.625) 

TEST=. FALSE, 

K=INX(KR»1) 

IZR=INDS(K) 

J=INX(KP*1) 

L«INDS( J) 

IF(L.NE.l) GO TO 70 
IS=J 

IF(KR.EQ.l) GO TO 6 

L=KR-1 

DO ^ I-=i,L 

A IS=IS+4*INDS( I+l) 

6 L=INDS(XR+1) 

AMR=DS(K+1) 

IF(KR.£Q,KP) 60 TO 40 
AMP=DS( J+1) 

ALPH=DS( J+4) 

RM=1.0/AMR+1.0/AMP 
rM=(T*RM)** 0.3 
LS=INX(i^R,4)-l 
CU:s( IZR-1 ) 

DO 30 I=1»NL 

K.= ( ( I-1)*(2*L-I > )/Z 

C6=BP( I+IS+L) 

IFUZR.EQ,!) GO TO 16 
C3=DS{ I+LS) 

C3=F(C3) 

C6=C6*CU*SORT(C3) 

16 DO 30 J=I,NL 

IF(J.EO.I) GO TO 20 
C6P=BP( J+IS+L) 

IF(IZR.EQ.l) 60 TO 18 
C3=DS( J+LS) 

C3=F(C3) 

C6P=C6P*CU*SQRT ( C3 ) 

18 C6F>=ALPH*A8S(C6-C6P) 

W( J+K) =1,9E-3*C6P**0.4»RM 
GO TO 30 
20 W(J+K)=0,u 

30 CONTINUE 
GO TO 100 

40 RM=S(3R‘t(2.0*T/AMR) 

DO 60 1=1. NL 

K={ ( I-l)*(2^tL-l ) )/2 

C3=BP(I+IS) 

CU=BP( I+IS+L) 

DO 60 J=I»NL 
IF(J.E(3.I> GO TO SO 
C3P=C3+bP( J+IS) 
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CUP«CU+BP( j+lS+U i 

C3P*0,0969*C3P 

CUPg2.04E-4*C UP»RM 

W ( J+K » = AMAX I ( C~3P • CUP ) 
GO TO 60 

50 WJ J+jC)=0,0 

60 CONTINUE 
GO TO loo 

._.7 0 TESTS. TRU E.g_ 

100 RETURN 
END 
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SUBROUTINE EWTHl TE»KR*NL) 

INCLUDE SPEC 1»LIST _ » 

paramete'r NH~SQ=(NLm+NSQ‘)/F 
dimension W(NHSQ) 

equivalence ( IS>T'MPU) ) * (L»TMP( 2) ) » ( L1.»TMP{3.L) # t L 2 1 TM P f 4 n i t L3 « 

C TMP(5» ) *<2*TMP(6) ) *(CS*TmP(7) )*(RS»TMP{8) »i(KS»TMP{ 9) ) ;tS* 

C TMP(lu) )»(TD.TMP(lin»lTR,TMP(12))»(|C»TMP113ntlCl»TMPaA) )« 

C ... ( C2 *TMP (15 U » ( C3 , TMPJ 16 ) | ♦ ( C4 » TMP. ( ij ) ) > ( CIP » tWP (.18 ) J » LC3P » 

C ■ TMP(19)),(C4P,TMP(20))»(W(1)»TMP(21)) 

DATA Al*A2,A3fA4/0*029»l2.2»5.625>7. 6/ 

I$sO 

iF(KR.EQ.l) GO TO 6 
L=KR-1 
DO A I=1»L 

A IS=IS+A*INDS( I+l) 

6 L=INDS(KR+1) 

IS=IS+L 
L1=INX(KR»1) 

L2=INX(KR,A)-1 
L3=INX(KR»6)-1 
Z=INDS(L1)-1 
Z^Z*Z 

CS=DS( Ll+5) 

RS=DS(Ll+6) 

KS*INDS(Ll+7) 

S=RS+0, 5 
TD=TE**S 
TR=SQRT(T£) 

S=S/( l.o+S) 

RS=TE**RS 
DO 50 I=1*NL 
K=( ( I-l)*(2*L-I ) )/2 
C3=DS( I+L2) 

C3=C3*(C3+1.0) 

C3=C3/(C3*AA-A3) 

C3=SQRT(C3)*8P( I+IS) 

CA=BP(I+IS+L) 

I F (I NOS ( I+L3 ) »G T • 0 *0R* I * £Q«1) GO TO 10 
C 1=CS*DS ( I+LI+7 ) 

C2=A2*( Ai*Z/Cl) **S/TD 
C1=RS*C1*AMAX1( 1.0, C2) 

GO TO 12 
10 C1=0.0 

12 CONTINUE 

DO 50 J=I,NL 
IF(J.EQ.I) GO TO AO 
C3P=DS(J+L2) 

C3P=C3P*(C3P+1.0) 

C3P-C3P/ ( C3P*AA-A3 ) 

C3P*SQRT(C3P)*bP( J+IS) 

C3P=ABS(C3-C3P) 

CAP=ABS(CA-BP( J+IS+L) ) 

IF( INDS( J+L3) .GT.O.UR.J.cQ.l) GO TO 20 
C 1 P=CS*DS ( J+L 1+7 ) **KS 
C2=A2»(A1*Z/C1P)**S7TD 
CiP=RS*ClP*AMAXi( l.C,C2) 

GO TO 22 
20 ClP=U.u 
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22 CONTINUE 

CXP»CI+C1P 

i C2« CBRT 1TR*^P*C4P ) 

• j:4P»ii*i4*(2*c4P)**o.4/tR 
C4P«0o 0362*AMAX 1 ( C4P ♦ C2 ) 
C»P«AMAX1(C4P>C1P) 
C2»8.18*CBRT'«Z*C3P*C3 P)/T'r 
C3P*0, 0969*AMAX 1 ( C3P *C2 ) 

VM J+K)gAMAXl( C3P » C4P ) 

GO TO 50 
40 W(J+K)=0.0 

50 CONTINUE 
RETURN ■ 

END 
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-TFOR*Si SUB7 

-..subroutine LlNWIO(TE»T»AN£>AN»i(;R>NU 

INCLUDE. SPE.CltL.IST __ 

include SPECAtUiST 
LOGICAL TEST 

_LS=0 

IR=0 

•v 1F<KR*EQ*1) go TO 6 
' LaK R-1 
DO A I = 1-*L 
IS=IS+A*INDS(i+i) 

A 1R =1R-HNDS( I+l) 

6 L*rNOS(K.R+l) 

IS®IS+3*L 

. .JCALL EWTH(TE*KR.*NL) 

DO 2 ^ I=1»NL 
■ KsIFN(i*L) 

DO 2u J=I*NL 
20 W( J+K)*ANE*rtO{J+K) 

DO 30 ksl^NSP 

CALL NWTH(T»KR»K*NL>TEST) 

IF(TEST) GO TO 30 
DO 25 1=1»NL 
N=IFNU »L) 

DO 25 J=I*NL 

W( J+N)=W( J+N)+AN(M*WD( J+N) 

25 CONTINUE 
30 CONTINUE 

N=INX(KK.2)-1 
A=INDS(N-7)**2 
A=1*A15E3*SQRT(T/DS(N-6) ) *A 
DO 50 l=i.»NL 
K=IFN(I »L) 

B=BP( I+IS)+R( I+IR) 

C=i.O/DS( I+N)**2 
DO 50 J=I*NL 

W( J+K)=W( J+K)+B+3P( J+IS)+R(J+IR) 

50 WD( J+K)=A*(C-1,0/DS(J+N)*^^2) 

RETURN 

END 
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TFOR»SI SUb6 //-xi 

SUBROUTINE cRAT E IT E »KR » NL ) 

INCLUDE SPEClfLIST 
PARAMETER N2=NSG+NLM+1 
dimension RM{NLM»NPL) »STP{15) 

EQUIVALENCE { RM ( 1 * 1 ) »TMP { 1 ) ) , C STP { 1 ) ,TMP(N2 ) ) , (L.STPd ) ) » (Ll» 

C STP(2) ) »(L2»STP(3) ) ♦(L3,STP(4) )» (L4»STP«5) ) , (ANG»STP(6) ) , 

C <U*STP(7) ) , tV,STP(8) ) ,tVP»STP(9) ) ♦ ( GL »STP ( 10 ) ) »{GU»ST?(11) ) , 
C <A*STP(12) ) »(B,STP(13) )*(Y»STP(14) ) *(R*STP(15) ) 

DATA PC/15*79/ 

L=INDS(KR+1> 

L1=INX(KR.2)-1 

L2=INX(KR»5)-1 

L3=INX(KR»6)-1 

L4-INX(KR,8)-L-1 

ANG=INDS(L3+1) 

Z=INDS(Ll-7)^^^^2 
DO 60 I-1,NL 
U= Z/DS( I+Lll»*2 
GL=DS( I+L2) 

DO 50 J=I.NL 
IFU.EQ.I) GO TO 4C 
V=VP < 

VP= Z/DS( J+.Ll+l )**2 
GU=DS( J+L2) 

A=U-V 

Y=A*PC/TE 

R=GU^<EXP(-Y)/GL 

IF{I,EQ,1.AND,IABS(INDS(J+L3) ),EQ.2) GO TO 30 
6 B=DS( I+L*J+L4) 

IF(B*LE*0.0) GO TO 20 
B=A«^SQRT(A) 

A=U/A 

RM( I »J)=ERFGRY( Y,A)/B 

A=U-VP 

Y=A*PC/TE 

B=A*SQRT(A) 

A=U/A 

RM(.I»JirRM( I»J)>ERFGRY(Y»A)/B 
10 IF(I.£Q*1) RiM( I .J)=ANG*RM{ I»J) 

12 RMIJfl )=RM( I»J)/R 
GO TO 50 
20 B*(U~VP)/A 
A=U/A 

RM( I »J) = ( A/yj *SQRT { A/U ) *EREG { Y » A » B ) 

GO TO 10 

3L B~ERFGS(TE. J»KR) 

_JEJ.B«LT#Q.01 60 TO 6 
R'M(I*J)=B 
GO TO 12 

40 _ 

VP= Z/DS( J+Ll+1 >»*2 
50 CONTINUE 

yg Z/D SINL+LH-l ) »*2 ■ 

A=U-V 

y=a«pc/te 

B=A*BQRT<A) 

A=U/A' 

RM(I*NPtt=ERF6RY(Y»A)/B 

r.42:^ .. 
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IF( I .EO.D GO TO 54 
GO TO 60 

54 VP=ERFGSCTE»L+1»KR) 

IF(VP.GT.C.O) GO TO 56 
RM( I »r<PL)=ANG*RM( I »NPL) 

GO TO 60 
56 Y=U»PC/TE 

B=U*SQRT(U) 

RM( I »NPL)=AMG»( RM( I ».\PL) -ERFGRYJ Y » 1. C )/B)+VP 
6C COF^TlNUE 
RETURN 
END 
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Functions TVEF and AXEF give transverse or side and axial or end escape 
factors for inverted transitions, respectively. Least square curve fits to 
numei'ical calculations are used (see Appendix B) . l!he arguments are C - axial 
gain logarithm, R - cavity aspect ratio and P - ratio of Lorentz to Doppler 
line widths. To avoid possible overflow problems, values of gain are limited to 
Q 0 i 7 ^ain bounds in the functions. The function CSA produces the average side 
solid angle fraction for L _< 0, the average end solid angle fraction for L = 1, 
<r/L> j , for L = 2 and <r/L> , for L > 3. Function LPI gives <p> and WLP 

gives the line profile constant k for the input of the line width ratio. 

The subroutine ESCAPE calculates the emission factors E^^ for the transition 
j i. (Cavity effects are not involved.) The factors use the symbol PHI and 
reside in TEMPS, Required input is the excited level populations ANX for type 
KR, The optical depth or gain logarithm is denoted by G. An expression needed 
for the Jacobian matrix, the rate of change of the logarithm of the emission 
factor with respect to the logarithm of G, is placed as the symbol W in TEMPS, 

If the magnitude of the excess emission (Appendix B) is less than 10 , this 
expression equals the excess emission and is not needed; W then contains a multi- 
plicative factor that describes the power broadening. The broadening is simply 
this power factor times the Einstein A-value (noncavity effects only). For the 
larger magnitude excess emission, the power factor equals the excess emission 
factor times the upper level population density and divided by the inversion 
density. (For absorbing transitions, both the excess emission and inversion 
density are negative giving a positive power factor.) Note that the excess 
emission and inversion density approach zero at the same rate so that the power 
factor as contained in W can be determined from an expression that is a limit 
ratio. Arrays GP and LU in TEMPS are used to store data on inverted levels in 
the following manner: GP(l,i) are the positive gain logarithms in order of 

descending values (up to i=4) , GP(2,i) are the corresponding line width 
ratios, GP(3,i) the transition center profile values, LU(l,i) the lower level 
indices and LU(2,i) the upper level indices. The data corresponding to a cavity 
transition are used in the cavity rate equations. If a cavity transition is 
not one of the four maximum of gains greater than unity for a given atom or ion 
type, then data are not made available for the cavity and coupling of the cavity 
Intensity to the laser medium is neglected. (The diagonal elements of the 
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emission factor are used to hold the total gain of a level for transitions to 
all lower levels . An upper bound on any neglected cavity coupling thus may be 
psrflbTished. 'I 
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FORtSi SUB9 

. FUNCTION TVEF(C.R.P) 

DIi4ENSION A(7f3 ) .B(3) 

DATA A/6, 9621 3. -l*3732f 4. 3505 »-4.768 3, 2.4813 »-0* 59608, 0*0&2728» 

3. 1943 ,-7.6577 ,12. 775, -12. 875, 6. 6044, -1.5975,0. 14369, 3i840l» 
3. 7817,-15-. 395, 17. 138, -8. 7 419, 2. 0886, -0.18833/ 

include SPEC1,LIST 
INCLUDE SPEC2,LIST 
REAL LPI 

X=AMINK3.912,AL0G(R) ) 

DO 10 1=1,3 
B‘U)-0.u 
DO 10 J=7,l,-1 
Bt I)=X*B{ I ) 

10 B(I)=B(I)+A(J,I) 

D=CFi5T(3) 

ENTRY TV2(C»R,P) 

CE=C 

IF(C.GT.50.0) CE=50.0+30.0*(C-50.0)/ (C-20.0) 

W=l. 0/(1. 0+12. 5»(CE/R)»»4) 

TV£F=( 1.0-W)»WLP(P)*T\/Q(CE)+W»CE*D*LPI (P) 

REfuRN 

FUNCTION T/G(OUM) 

X-=0.1*DUM 

X=-B( 1)+B(2)»X+B(3)*X*X 
IF(X.GT.IC.O) X=10.0 
TVQ=EXP(X) 

RETURN 

END 
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FUNCTION AXEF(C*R*P) 

DIMENSION A(6,3) »B(2) tCA(3) ‘ • ' L. 

DATA A/-8.47794E-3,l,20i25E-l*5*0607 3E-3.-6,869'03E-3*8.67575E-^ 
C -3,25246E-5*-1,59631E-1»2.80494E-1.“1.21696E-1»2*19665E-2» 

C -1*79975E-'3»’5«5E-5,-X,42061E-2»1.06835E-1»“4*34725£-2j! 

C 6«91377E-3»-4*9476E-4,1,33524E-5/B/0.0337633»-0*0080375/ 

INCLUDE SPECltUIST 
INCLUDE SPEC2>LIST 
REAL LPI 
CE=C 

IF(C.GT*5G,0) CE=5C,0+30*0*(C-50.Q1/ {C-20.0) 

DO 10 1=1.3 
CAU )-0.0 
DO lU J=6»l»-1 
CA(I )=CE*CA{I) 

IvJ CA( n=CA( I )+A( j.n 

IF(CE*LE.4.0) GO TO 20 
G=D( 1)+B{2)*CE 
IF(CE.GE,6,0) GO TO 19 
W=2*0**CE/16,0 

CA{2)=CA(2)+(1.0-W)*(CA(2)-G)/3.0 ’ 

•GO TO 20 

19 CA12)=G 

20 G=AL0G(R) 

IF (CE.LE.1.0) go to 40 

W={CA{ 1)+CA(2)»^G+CA{3)*G*G}/(CE*CE) 

30 G=W^^EXP(CE)*SQRT(CE)/{R^^R) 

GO TO 50 

40 W= 0,1 10674+0. 01 9589*G+0.055589*G*G 
GC TO 30 

5C W=1.0/ ( 1,0+50. 0*CE) 

AXEF={ l.C-W)»G*WLP(P)+W»CE*CFST (4)*LPI (P) 

RETURN 

END 
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-TFORiSl' SUBil 

FUNCTION- CSA(R»L) 

___ _DJMENSlON A14.2) 

DAtA" A/0,42-44.-0*172#0.09*-0. 025*0. 66, -1.09*2. 3, -1«6/ 

IF(L.GE.3) 60 TO 30 

M^MAX0(1,U 

B=*0.0 

-DO IC I=4*l*”l 
10 _B={B+A{L»M) )/R 

IFiL.GT.O) GO TO 20 
CSA=*1.0-B 
RETURN 
20 CSA=B 
RETURN 

30 CSA= (0.9428+ALOG (R) ) / ( 8 ,0*R*R) + 0. 0132/R^^*4 
RETURN 
END 

-FOR* Is SUBUB 

FUNCTION LPI(X) 

REAL LPI 

LPI=(0.707+1.25*X)/( 1.0+2. 5*X) 

RETURN 

END 

-F0R*IN SUBllC 

FUNCTION WLP(X) 

WLP=(1.0+Q,7304*X+0.5811*X*X)/( 1.0+1.2946^‘X + 1.0299*X*X) 

RETURN 

END 
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-TFORiSr ’ SUB12 
-i • SUBROUTINE ESCAPE(ANX*KR*NL*AN) 
sPEiifkisT 
INCLUDE SPEC2*LIST 
INCLUDE SPECA»LIST 
dimension ANX(NLM) 

REAL LPI 

FIX) = (1*0+0*936A4*X+0*3299*X«-X) /{ 1.0 + 2.06A82*X+1.65979«X*X+ 

C.^ ^58473tX*X*XJ 

DEP(X) = (srGN{0.25»X)-0.75) 

L=INDS(KR+1) 

RLtL— . 

RL=RL+0.5 
L1=INX(KR»5 )-l 
L2=INX(KR»8)-1-L 
KD=(L*(L+l)')/2 
DO 10 K=l»4 
6P(1»K)=0.0 
GP(2>K)=0,C 
GP(3»K)=0,0 
LU(1»K)=0 
10 LU(2»K)=0 

DO 100 K=KD»1>-1 
P = 2*K 
P=RL»RL-P 

I=RL+1.0-0,5*(SQRT (P)+SQRT{P+2,0) ) 

J=K-IFN( I »L) 

IF{J.6T.NL*0R.J.EQ.I) go to 90 
P=DS( J+Ll )/D5U+Ll ) 

RS=1.49736E5*ABS(DS(J+L*I+L2) )*CL/WD(K) 

G=RS*( ANX{ J)-P»ANX(I ) ) 

P=W(K)/WD(K) 

PF=F(P) 

IF(G.GE.O.O) GO TO 50 
PA=-G/(2.0*CAR) 

PF=PF*LPI (P)*(CFST(3)+CFST(4) ) 

IF( (-G»PF).LT.l*0E-3) GO TO 20 
G=-G^tPF 

PHI (K)=0.84*SQRT { P/PA) 

PA=C.9/(PA*SQRT(AL0G(AMAX1(PA»2.72) ) ) ) 

IF(PA.GT.PHI (K) ) PHI(K)=-PA 
RS=ABS(PHI (K) ) 

PA=5.0*RS 

C ARRAY W IS USED IN JACOBIAN MATRIX 

WnO = ( 1*0+PA-PA*PA)">EXP{-PA)*RS»DEP{PHI (K) )-6*EXP(-G) 

PA=R5*( l.G+PA)*EXP(-PA)+EXP(-G) 

W(K)=W(K)/PA 
PHI(K)=PA 
GO TO 3u 

20 PHI {K)=1.0+G»PF 

WIK)=RS*PF*ANX( J) 

GO TO ICO 

30 IF(ABS(PHI (K)-l.O) .LT.l.OE-3) PHI ( :<) = 1 .0+SI GN ( 1 . OE-3 ♦ Pril { K ) -1 .0 ) 
GO TO 100 
50 G=G*PF 

RS=RS»PF 
DO 60 Jl=l»4 

IF(G.LE.GP( 1»J1 ) ) GO TO 60 
DO 54 J2=4,l*-1 
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IF(J2.EQ.J11 GO TO 56 
GP(l.J2)=GPll,J2-l) 

GP(2tJ2)=GP(2,j2-l) 

GP{3tJ2)=GP(3,J2-l) 

LU(l»J2)^LUa,J2-l) 

LU(2.J2)=LU12,J2-1) 

54 CONTINUE 
56 GP(1.J1)=G 

GP(2»Jl)=P 

GP{3»J1)=0.5642*PF/WD(K) 

LUdf J1) = I 
LU(2»Jl)=J 
GO TO 62 
60 CONTINUE 
62 J2 = l- 

DO 70 Jl=i,NCAV 
IFUCAVUtJD.NE.KR) GO TO 70 
IF( ICAV(2»J1) .NE. I ) GO TO 70 
IF( ICAV(3»J1) .NE.J) GO TO 70 
J2=0 

70 CONTINUE 
PA=CFST{3) 

■ IF(J2.GT.O) PA=PA+CFST(4) 

PF=' PA*LPI(P) 

IF( (G»PF) .LT,l,0E-3) GO TO 20 
PHnK)=TV2(G»CAR»P)’ 

IF(J2*LE,0) GO TO 80 
PHI (K)=PHi (K)+AXEF{G,CAR»P) 
RS=TV2(1.01*G»CAR»P)+AXEF( 1,01*G»CAR»P) 
GO TO 85 

80 RS=TV2{1.01#G,CAR,P) 

85 PA=PHi(K) 

_RHI (K) = 1.0+PA 

W ( K ) = 1 0 0 , 6*P A* A LOG ( RS/ P A ) >PH I ( K ) 
J1=J+IFN( J»L) 

RS=AM1N1(G»88,‘0) 

- PHnJl)=PHI{Jl)+EXPlRS) 
go’ to 30 

90 PHI (K)=0,0 

100 CONTINUE 
RETURN 
END, 


- 55 - 



77--11 


The subroutine EQUIL is used to establish an equilibrium initial excited 
level distribution. Let 

- L . 

n 

S Z S 1 

then the total population of type s,z is ^ • 

If 


P 


s ,z 

* s 


1 , 


p 


S 3Z 


then 


z 

s 



y=z+l 




3 

Specie conservation is given by N « N- Q where 

S X s 

S z -z+1 

(N^) " a®-" . 

z 

Charge conservation is 

S “s'*. ■ L ) • 

s ' S 

Function Q Is a polynomial in N^. . An initial estimate of is obtained by 
computing the next stage ionization fraction of each type ordered with Increasing 
3j^. A fraction ^0. 86 is considered equivalent to unity; when the fraction is 
less, the calculations are halted. The N estimate is used to calculate 
95,11 Q /dZn N and hence a new N , the new and old estimate are averaged and used 
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to calculate a new N , etc., until convergence is obtained. Once N is known, 
populations are readily calculated. The symbol ANG is used for ground state 
populations in the subroutine rather than AN. STOR, AN and ISZ are scratch 
variables used for ordering? other vector quantities are self explanatory or 
conventional . 

The subroutine QSTAT is used to establish a quasistationary initial distri- 
bution. The form of this distribution is 0 if and only if z z^(s) and 

z = 2 ^+l, n=* 1. That is, only one ionization stage of each specie is con- 

side red in detail • If 


s , Zj+1 

=Vs • 


then 


(1 


S,z 

- a )N = / , K 

3 3 ^ 

n 


and 


"a ■ S - I + “s>"a • 


Inputs through namelist NAME14 are z^ as IZ and P to determine a^. If.R is 
negative, 0 < |p| f 1. then ct^ is set equal to |P| for all s. If P is positive. 


a 

s 


(i + f \ ' • 


P=1 corresponds to equilibrium and P*'! for a recombining plasma. Levels for 
ji>2 are assumed quasistationary. For positive P, N^ must be obtained by an 
iterative procedure. If z^ > z^ or z^ £ 0, then specie s is considered fully 
ionized into the z^+l ground level. Populations are determined by iteration 
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using subroutine SRJM and starting with essentially null quasistationary levels 
“8 

(value of 10 )• Iteration Is halted when the, change of ground level populations 

of are sufficiently small. 

EQUIL and QSTAT place ground state populations (AN or ANG and ANP) into 
MSTOR common for use by the program at integration start. 
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SUBROUTINE EQUIL 

INCLUDE SPECl.LlST 

include SPEC2.LIST 

INCLUDE SPEC7*LIST 

PARAMETER N13S0«l+2*( NSP+NCM> 

PARAMETER NIT«10*N13Sl«NSP+l4tN13S2=N13Sl+NSP*N13S3*N13S2+NSP» 

C N13S4»N13S3+NSP *N1 3S5»N13S4 +NSP 

0 1 mens ion bet a ( NSP ) *S I GMA ( NSP ) , P ( NSP ) * AN { NSP ) * I SZ ( NSP ) * STOR ( NSP ) 
C AN6(NSP)fANP(NCM) 

EQUIVALENCE ( ANR ( 1 ) vSTRMTX (!))»( ANPM ( 1 ) tSTRMTX ( 1 > ) H ANG ( 1 ) • 

C STRMTXI 1 ) ) * ( ANP(l) »STRMTX<N13SO) ) 

EQUIVALENCE (LtTMPd) ) * (LN*TMP( 2 ) ) » ( LG*TMP( 3 J ) » ( MC,TMP (4)) , 

C (A*TMP(5) )♦ (B»TMP(6) ) »(CtTMP(7) ) » (D*TMP( 8) ) » (E*TMP(9) ) > 

C <RI#TMP(XO) )»(GNS*TMP(lin* (Q*TMP(12n »(S>TMP(13)) * 

C (BETA(l)»TMP(14n »(SIGMA(1) »TMP(N13Sl) ) * (P( 1 ) *TMP(N13S2) » ♦ 

C (ANU)*TMP(N13S3) ) t (1S2(1),TMP(N13S4)) f (ST0R<1>» 

C TMP(NI3S5) J »(CGMEfTMP(l))»(DGME*TMP<2))» (GLTEtTMPO)) 

4 FORMAT <1H0*2X30HSUM ON 2 IN SUBR EQUIL - ERROR) 

5 FORMAT (1HO»2X16HSLOPE FUNCTION =1PE10.3»2X10HREL ERROR=1PE10. 3» 

C 2XI2»1X13HSTEP IN EQUIL) 

6 FORMAT ( IH0»2X9HF INAL NE»lPEiO-3.2X6HERROR=lPE10*3) 

DATA ERROR/1. OE-3/ 

Bs2.0708E-7/(TE*SQRT(TE) ) 

C=15.79/TE 
DO 10 K«1,NSP 
LN=1NX<K*2)-1 
LG*:INX(K*5)-1 
L*INX(K»1) 

GNS=DS(L+2) 

L=INDS(L) 

RI=L 

L»INDS(K+1) 

SI6MA(ia»0.0 
DO 10 J«1*L 
A=RI/DS(J+UN) 

A»B*DS ( J+LG > *EXP ( C*A* A ) /GNS 
IFIJ.NE.l) GO TO 8 
BETA<K)«A 

8 S16MA(K)«SIGMA(K)+A 
id CONTINUE 

IF(NSP.EQ.l) 60 TO 17 
LN«1 
MC»1 

DO 16 K*2*NSP 
L=INX(K*1) 

L“IN0S(L) 

IF(L.EQ.l) 60 TO 12 
GO TO 15 
12 RI-1.0 

DO 14 J=MC»LN#-1 
PU)=RI 

14 RI=RI*BETA( J) 

LN*K 

15 MC=MC+1 
IF(MC.EQ.NSP) GO TO 12 

16 CONTINUE 
GO TO 18 

17 P(l)*1.0 
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18 ISZ<1)«1 

AN(l)*ANca) 

ST0Ra)*8ETA(l) 

IF(NSP.EQ.l) GO TO 26 
MC»1 

DO 24 K«2*NSP 
L*INX(Kfl) 

L«INDS(L) 

IF(L*EQ*1) MC*MC+1 
R.I*ANC(MC) 

AN(K)=RI 
STOR(IO*BETA(IC) 

ISZ(K)«K 
DO 20 J»|C»2#-1 
IF(BETA<K).GT*ST0R{J-1) ) GO TO 24 
ISZ(J)»I$Z(J-1) 

ISZ( 

AN(J)*AN<J-1) 

AN(J-1)»*RI 

ST0R(J)«ST0R(J-1) 

ST0R(J-1)*BETA<K) 

20 CONTINUE 
24 CONTINUE 
26 ’ GNSaO.O 

DO 40 J*1*NSP 
RI*AN( J)*STOR(J) 

B*GNS*STORt J)+1*0 
0*s4.0*RI/(B»B) + l*0 

IF(B*LT«1*09,AN0.D*LT.1.30) GO TO 36 
A* B*(SQRT(D)-1.0)/<2#0*RI) 
ANE=»A*AN( J)+GNS 
60 TO 42 

36 6NS*GNS+AN( J) 

40 CONTINUE 
42 NSTP=1 
44 C«0.0 

Dal*0 
LN=1 

DO 60 J#1*NC 
GNS»ANC(J) 

IF(J«EQ*NC) GO TO 48 
LG=IC( J+D-l 
GO TO 50 
48 LG*NSP 
50 L»INX(LG*1) 

LsINDS(L) 

RI»L 

C»C+RI*GNS 

MC«0 

RI»0.0 

0 = 1*0 

DO 54 K=LN*L6*1 
MC*MC+1 

B=SIGM A < K ) *P ( K ) *ANE** « L-MC ) 

Q=Q+ANE*B 
A*L-MC+1 
54 RI=RI+A*B 

D=D+RI*GNS/Q 
IFCLG.EQ.NSP) go TO 62 
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LN»L6+1 

IF(L.NE*MC) WRITE(6*4) 

60 CONTINUE 
62 GNS=C/0 

E*(GNS-ANE)/GNS 
IFINSTP.EQ.l) GO TO 64 
IF(ABS(E),LT,1.0E~8I GO TO 64 
S»D+2*0*I0-S)/E 

IF{NSTP.6E*NIT*OR*S.LE«0*0) WRITE{6#5) S*E»NSTP 
64 S*0 

ANE»0*5*(GNS+ANE» 

IF(A8S(E)*LT,ERR0R*0R,NSTP*GE*Nm GO TO 70 

NSTP*NSTP+1 

GO TO 44 

70 LN«1 

IFINSTP.LT.NIT) GO TO 71 
IF(ABS<E),LT.0.01) 60 TO 71 
NSTP*NSTP+1 

IF(NSTP«LT*3*NIT) 60 TO 44 
WRITE(6*6) ANEtE 

71 CONTINUE 

DO 80 J*1»NC 
6NS»ANC(J) 

IFIJ.EQ.NC) GO TO 72 
LG=IC( J+11-1 
GO TO 74 

72 UG«NSP 

74 L»INX(L6#1) 

L«IN0S(L) 

MC*0 

Qal*0 

DO 76 K*LN*L6*l 
MC=MC+1 

76 Q=Q+S 1 6MA ( lU *P ( K ) *ANE** { L-MC+1 ) 

MC«0 

DO 78 K*LN*LG»1 
MC=MC+1 

78 BET A ( K ) «6NS*P ( K ) *ANE** ( L-MC > /Q 

IF(LG.EQ.NSP) GO TO 82 
LN»LG+1 
80 CONTINUE 
82 MC=0 

DO fO K»1»NSP 
L*INDS{K+IJ 
GNSaBETAIlO 
DO 90 J«1*L 
MC«MC+1 

90 ANB{MC)=GNS 

CALL CONSRV ( 6T • ANP »G1 » G2 » G3 »G4> ANG • 0 ) 

E«<ANE-GT>/GT 

ANEsGT 

WRITE(6>6) ANEiE 

CGME*0.0 

DGME«0.0 

CALL ELCN(GTt0«0>0*0»0«0»ANG»ANP) 

GLTE=GT 

RETURN 

END 
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SUBROUTINE QSTAT 
J,L. =. IMCUUPE SPE.Ci#LlST_ 

INCLUDE. SPEC2. LI ST 
INCLUDE SPEC? »L 1ST 

DIMENSION NBRS(NSP) >ALPH(NCM)»BETG(NCM) »IS(NCM) »KS ( NCM ), » IZSJ NCMl .» . 
C’ ■ i ANP(NCM)»AN('NSP),R3(NLM) ,LC(4)*D(6) ♦RilNLM.NPL) .R2(NLM*NPL) • 

C R4(2»4)»R5{4) »R6{3»4)fI2(NCM) 

LOGICAL L GP»LERR __ 

PARAMETER N1=2*NSP+NCM+1 ♦N2=Nl+5*NCM+6 *N3=N2+2* ( NLM+NSQ ) ♦ 

C- - N14S1^*NSP+1»N14S2=N14S1+NSP,N14S3=N1+NCM»NI4S4*N14S3+NCM» 

C N14S5=N14S4+NCM,N14S6=Nl4S5+NCM»Nl-4S7=JS14S.6+NCM»Nl4S8=N2fNLM+ 

C NSQ*N14S9=N3+NLM,NIT=10 

EQUIVALENCE (CGME»TMP(1 ) ) * (DGME*TMP(2 ) ) ♦ ( GLTE» TMP ( 3 ) ) * ( AN( 1) » 

C STRMTX(l) ) » (NBRSd) »STRMTX( N14.S.1 ) ) liALPb ( I.) tSTRMTX ( N14.S2J.) » . 

C ..{BETGIl') jSTRMTXC'nI) ) » (ANPd ) »STRMTX ( N14S3 > ) » { I2S( 1) » 

C STRMTX(N14S4) ) » ( IS( 1) ,STRMTX I N1455 ) ) t (KSIl ) »STRMTX ( N14S6 ) ) * 

C (D(l) >STRMTX(N14S7) ) ♦{R1(1*D >STRMTX(N2) ) » lR2(lfl) » 

C STRMTXIN14S8) ) >{'r3(1) ,STRMTX(N3) ) ♦ ( R4 1 1 , 1 ) >STRMTX( N14S9 U * 

C (ANRI 1) »STRMTX(1) ) 

NAMELIST/NAME14/EBROR,P,IZ 

1 FORMAT! 1HQ,3X18HQSTAT ERROR, TlTPE 12) 

2 FORMAT! 1H0,3X14HQSTAT NE DI FF=1PE10. 3 ) 

3 FORMAT! 1H0,3X14HQSTAT NG DI FF=1PE10. 3 ) 

4 FORMAT !1H0,3X15HQSTAT GIVES NE=1PE10.3) 

READ(5»NAME14) 

LI0=0 

LGP=. FALSE. 

DO 10 K=1»ND 
BT!K)=0.0 

10 AN8(K)=1.0E-8 
DO 11 K=1,NSP 
AN(K)=0.0 
NBRS(K)=NBAR!K) 

11 NBAR(K)=2 
DO 12 K=l,4 

12 ANPH!K)=0.0 
C=2.0708E-7/!TE*SQRT(TE) ) 

IF(P.LT.O.O) GO TO 50 
DN=0.0 

A=15.79/TE 
DO 20 J=1,NC 
AiNP! J) =0,0 
IS! J)=0 
KS(J)=0 
K1=IC! J) 

IF(J.LT.NC) K2=IC(J+1)-1 
IFU.EQ.NC) K2=NSP 
MK=K2-K1+1 

IFIIZ! J) .GT.MK) 60 TO 18 
IF! IZ! J) .LE.O) GO TO 19 
DO 15 K=K1,K2 
MK=K-K1+1 

IFIMK.NE.IZU) ) GO TO 15 
Z = MK 

MK=INX(K,l)+2 

GD=DS(M.<) 

MK=INX(K»5) 

GD=DS(MK) /6D 
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Z=Z/DS(MK) 

BETGt J ) =C*GD*EXP (A*Z«Z) 

B=IZ( J)-l 
60 TO 

15 CONTINUE 

18 IZ{J)=-IZ(J) 

19 BETGIJ)=0.0 
B=MK-1 

20 ON=DN+B»ANC(J) 

IT=C 

ANE=DN 
2A Y=0.0 

DO 30 J=*1,NC • 

30 Y=Y+ANC ( J) / { 1.0+P*ANE*BETGt J) ) 

A=DN+Y 
IT=IT+1 
B=( ANE-A)/A 
ANE=(ANE+A) /2.0 

IF(ABS(B)*LT. ERROR. OR.IT.GT. NIT ) GO TO 32 
GO TO 24 

32 IF( IT.GT.NIT) GO TO 36 

33 00 34’'j = l»NC 

34 ALPH{ J) =1.0/( 1.0+P*ANE*BETG( J) ) 

GO TO 60 

36 WRITE(6»2) B 

IF(ABS(B).LE.0.01) GO TO 33 
IF(IT.LT.2*NIT) 60 TO 24 
IF(ABS{8).GT.0.05) go TO 41 
GO TO 33 

40 WRITE(6»1) LIO 
STOP 

41 LI0=1 

60 TO 40 

42 LI0=2 

GO TO 40 

44 LI0=4 

GO TO 40 

45 „LI0s5 

'go to 40 

50 A=0.0 

IF{ABS(P+0,5) ,GE.0.5) 60 TO 42 
DC ”56 J = lfNC 
ANP( J)=G,0 

KS(J)=0 

ALPH(J}=-P 

iFjO«LT.NC.) K2=IC(J+1)-1 
IF(J.EQ.NC) K2=NSP 
IT=K2+1-IC( J) 

IFUZt J).GT*IT) IZ(J)=-IZ(J) 

IF( iZ( J) ,6T.0‘) GO TO 52 
B=IT 

GO TO 56 
.52 B=IZU)-1 
8 * 8 - 1 =: 

^6 A? A+B *ANC(J) 

ANE=A 

60 DO 61 J=ltNC 
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61 ALPH( J)=ALPH( J) *A.NC(J) 

MK=C 
QM=C *C 

DO 6 A J=1,.'K 
iCl=IC( J) 

ircJ.LT.NC) K2=IC(J^lj-l 
IF(J.EQ.NC) K.2=NSP 
IZS( J)=K2-K1+1 
B=IZS( J) 

QM=QM+B*ANC( J) 

DO 63 K=K1»K2 
IT=K-K1+1 

IF(IT.NE.IZ( J) ) 60 TO 63 
KSt J)=K 
IT=MK+1 
ISC J)=IT 

ANIK)=ANC( J)-ALPH( J} 

BTUT)=BET6(J) 

ANB(IT)=AN{K)/( ANE*BETG( J) ) 

IF(K.EG.K2) GO TO 62 
AN(K+1)=ALPH{ J) 

A=15.79/TE 
Z=IZ( J)+l 
IT=INXC<+1,2) 

Z=Z/DSl IT) 

IT=INX(K+1,5) 

GD=DS( I T ) 
lT = INX(iC+l,l)+2 
6D=GD/DSC IT) 

B=C*GD*EXP( A^i-Z«-Z) 

IT=IS( J)+INDS(K+1) 

BTC IT)=B 

ANBC IT )=ALPH( J)/ C ANE#B) 

GO TO 63 

62 ANPCJ)=ALPHC J) 

63 MK=MK+INDSCiC+1) 

64 CONTINUE 
IT=t^ 

70 DN=O.U 
A=0,Cj 
B = 0.u 
C-O.C; 

DO 80 J=1»NC 
IF( IZC J),LE.U) go to 78 
ANHS=0.0 
K=KSC J) 

MK=INDSCK+1)+IS( J) 

BTPL=BT(MK) 

IFUZSCJ) .NE.IZC J) ) GO TO 71 ■ 

ANHS=ALPHC J) 

BTPL=l.o/ANt 

71 MK=ISCJ)-1 
ZS=iZSC J) 

Z=IZt J) 

GD=ANBCMK4-1 J 

anep=ane 

CALL SNJMCRi»R2»R3»AI »AR»ARR*GRI ♦QRR »K»MKtAN»ANEP »QM»ANPH,ANHS» 

: BTPL*u,O»0.0*C.0»ZS*Z»DC 1) , LC »R4 » LGP » J » 0 » D C 2 )»D(3)jp(A)»D(5)» 

; DC6) ♦LERR»Ki»R3»R6) 
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IF(LERR) GO TO 44 
ZS=ZS+l-0-Z 

A= A+ZS* t ANE*A I - ALPH ( J ) * AR ) 

B=B+ZS*ALFri I J ) * ARk 
C=C+ANt*( ALPH{ J )*QRR-ANE*uRI ) 
mk;=is( J) 

ANHS=0*o 

Kl=NSTAR(K.)-2 
IF(Kl*LT.l) GO TO 73 
DO 72 K2 = l,»a 

72 ANHS=ANHS+AN£*BT ( MK.+K2 ) *ANB ( MK+K2 ) 

73 ANHS=ANC( J)-ALPH( J)-ANHS 
IF(ANhS.LT.O.O) ANhSsC.O 
ANHS=ANHS/(AiNE*uT(MK.) ) 
ZS=<GD+ANH6)/2.0 
IFtZS.EQ.O.u) GO TO 75 

Z =ABS{60-ANriS) /ZS 
IF(Z.GT*DN) DN=Z 
75 ANB(MK)=Zb 

AN(IC)=ANE*ZS*BT (MK.) 

78 IF(J.LT*NC) GO TO 80 
iT=IT+l 

IFCDN.LT. ERR0R.0R.it. GT. NIT) GO TO 82 
8U CONTINUE 
GO TO 70 

82 IF( IT.GT.NIT) GO TO 84 
GO TO 86 

84 WRITE(6»3) ON 

IF(DN.GT.0.05) GO TO 45 
86 WRITE(6»4) ANE 
DO 90 K=1»N6P 
90 NBAR(K)=NBRs(K) 

CALL ELEN(DN,O.O.C»A,AN,ANP) 

CGME=A 
DGME=A+B 
GLTE=ON 
•RETURN 
• END 
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Subroutine CONSRV calculates the electron density and densities of maximum 
ionization given the electron temperature and excited level normal populations. 
Densities of ground levels are also output. Symbol EPD is used for the electron 
density, EP is the equilibrium bound states factor e, DI = is proportional 

to the lowering of the ionization potential, DEL = A/N^ - B (see conservation 
discussion in main text) is a measure of the correction term for finite e in N^, 
QM = (1 + D)Ng is needed by subroutine SRJM and IBT is set £0 if array BT of 
Saha factors is to be calculated and >0 if not (e.g., during a second pass at a 
particular electron temperature) , 



n n 
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-TFORjSI SUB15 

, SUBROUTINE CONSRV ( EPD » ANP * EP»DI i DEL* QM» AN , I BT ) 

INCLUDE SPEC1*LIST 
include SPEC2,LIST 

PARAMETER Ni5Si=NCM+9»iNl5S2=Nl5S]+NCM>Nl5S3 = N15S2+NCM»Nl5S4= 

C N15S3+NCM 

DIMENSION ANP(NCM) »SA ( NCM ) »SB( NCM ) »SC { NCM ) ,SD(NCM) »ZC ( NCM ) *AN ( NSP ) 
equivalence (MS*TMP(1) J »(MR*TMP(2)). (GP»TMP{3) ) * (ZC . TMP ( N 15SA ) ) » 
(A»TMP(5) ) ♦ <8»TMP(6 ) ) * 1<L*TMP (7) ) ♦ ( D *TMP ( 8 ) ) » ( SA ( 1 } *TMP(9) ) » 
(SB(1) *TMP(N1531) ) » (3C(i) »TMP(N15S2) ) »(SD( 1) *TMP(N15S3) ) 

DO 2 J=1*NCM 
• SA(J)=0.0 

SB( J)=U.O 
SC( J)=U.O 

2 SD(J)=u.O 
MS=0 
MR=0 

iFUBT.GT.O) GO TO 3 
A=15.79/TE 

C=2*07ObE-7/(TE*iQRT(TE) ) 

3 DO 10 K=1*NSP 
Li=INX(K*2)-l 
L2=.INX(K*5)-1 
L=INX( K»1 > 

GP=DS(L+2) 

L=IND3(L) 

IF(L.EQ.I) MS=MS+1 
ZC(1)=L 
L=1ND3(K+1) 

DO U; J=1*L 
MR=MR+1 

IF(J.GE.NSTAR(>0 ) GO TO .10 
IF(IBT.GT.O) GO TO 4 
B=ZC(1)/D3( J+Ll) 

B=C*DS ( J+L2 ) *EXP I B ) /GP 

BT(MR)=B 
GO TO 6 

4 B=BT(MRJ 
6 B=B*AN6lMR) 

D=ZC(1)-1*0 
SA(M3)=3A(M3i+B 
SC1MS)=3CCMS)+d*D 
IF(J.NE.l) GO TO 1C 
AN(K)=B 
D=D*3QRT(D) 

SB(M3)=3B(M3)+B*D 
SD(MS) *SD( MS) +0*0*1 ZC‘( 1 J-2*0) 

10 continue 

C=0.0 
D=1.0 

DO 2u J*1*NC ^ 

IF(J.EQ.NC) 60 TO 12 
K=1C{ J+D-l 
GO TO 14 
12 K=N3P 
14 L=INX(K»1) 

L=INDS(U 
ZC(J) = L- 

C=C+ZC(J)*ANC( J) 
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D=D+ZC( J)*SA( J)-SC(J) 

2 ^ CONTINUE 
EPO=C/D 
EP=IE+TE*TE/T 
EP=EPD/EPiHt3 
EP=S(JRT(EP) 

EP=1.3E-2*SQRT{£P) 

A=0,0 

B=-O..U 

DO 30 J=1>NC 
GP=ZC( J)*SQRT(ZCU) ) 

ANP( J)=l,o+EP*GP 
A-A+GP*ANC { J ) /ANP ( J ) 

B=B+SD.{ J ) -ZC { J ) *SB { J ) +GP* ( SA ( J ) +EP*S 3 ( J ) ) /ANP ( J ) 
30 CONTINUE 

EPD=(C-EP*A)/(D-EP*B) 

DEL=A/EPD-b 

QM=0*EPD 

DI=4.8bE-A*SQRT (EPO/TE+EPD/T) 

DO 40 J=1»NC 

40 ANP( J) = (ANC(J)-EPD*(SAU)+EP*SB( J)') )/ANP{ J) 

DO 50 K=1»NSP 
5u ' AN(K)=EPD*AN{K) 

RETURN 
END ' 
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Subroutine SRJM is used to compute rates and the Jacobian for a particular 
atom or ion. The argument list is as follows, 


RL 


JCB 


F 

AI 

AR 

ARR 

QCRI 

QCRR 

KS 

MS 

AN 

ANEP 

QM 

NPH 

ANHS 

BETAPL 

GLTE 

CGME 

DGME 

ZS 

ZC 

ANPL 

IC 


rate matrix for the reduced set of equations. The vector G^ used 
to couple to the next ionization stage is the set of elements 
in RL with second index NPL. 

Jacobian matrix for the reduced equations. The element set with 
second index NPL is the vector 

((z^ + I - z)\) ORg/9V . 

used for cross terms in the overall Jacobian. 

■ + 

the vector 3R /3N .plus (1 + D) 3R„/3N if z = z , also used for 
cross terms 

ionization coefficient 

recombination coefficient (collisional and radiative) 

radiative recombination coefficient 

coefficient of energy absorption divided by kT^ 

recombination heating, coefficient divided by kT^ 

atom or ion type index (input) 

location of ground level in ANB, BT vectors 

ground level density vector (input) 

predicted value of (input) 

(1 +• D)Ng (input, may be obtained from CONSRV) 
photon density vector (input) 

Nj (input) 

(input) 

(dZdt)ilnTg (input) 

predicted (input) 

N 3R /3N (input) 
e e e 

maximum specie charge number (input) 
parent ion charge number (input) 

(output) 

if LC(j) = k > 0, then index k is used in the GP and LU data 
arrays in TEMPS common for cavity equation j 
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CVP 


LGP 

NCORE 

NEXP 

KPEN 

RPEN 

AHMET 

ANMEQ 

EMET 

LERR 

IGMX 


FMD 

CVOUT 


cavity rate matrix: CVP(1,J) is the induced transition rate 

for cavity j and CVP(2,j) is .the total spontaneous emission rate. 
The equation is 

I^Npha) = CVP(l,j) + rig CVP(2,j) . 

logical constant set to true for consideration of the Penning 
effect 

specie index (input) 

index for excited level of Penning effect (set negative for donor) 
Penning rate constant (set negative for two electron excitation) 
net Penning rate (output) 

metastable density at first input (to receptor atom) 

' equilibrium metastable density at first input 
metastable excitation energy divided by kT^ at first input 
logical constant output, set true in case of a’ matrix singularity 
that prevents an inversion (quasistatlonary levels) 
location of maximum positive noncavity gain in TEMPS data arrays 
GP and LU. Set null if there are no inversions or if all inver- 
sions are coupled to the cavity. 

vector containing the mode spacing parameters P^Av 
cavity output array: CVOUT(l,j) is the gain logarithm and 

CV0UT(2,j) is the transition width ratio for cavity j 


The function Ej(x) is used in the calculation of radiative recombination; 
t is approximated by the define procedure FX, At the start of the subroutine, 
he collisional rate coefficient matrix is placed in JCB, Values of normal den- 
ity for the atom or ion ate placed in ANR. Corresponding Saha factors are 
laced in BETA. Array AM and SRMR are used to hold derivative forms of the rate 
atrix with respect to populations. 

The internal subroutine CAVITY calculates the cavity rate- coefficient 
C = s«/t and its logarithmic derivative with respect to the gain logaritm, RCP . 
he internal subroutine INVERS calculates the inverse of an upper quadrant of 
atrix AMS (corresponding to the quasistationary level indices) and multiplies 
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the inverse into the array BVS composed of other elements of AMS and possibly 
vector F. These operations serve in the formation of reduced rate and Jacobian 
matrices . 

The vector 3N /alj is proportional to the vector of Saha factors; it is • 

6 XI 

used in calculating cross terms involving different atom or ion types. The 
'reduced’ form of the Saha factors is output through STORE common in symbol ANR. 
Also at output, GP(3,i) contains the relevant transition escape factor for 
i > IGMX > 0, and for cavity transitions contains the upper level to ground 
escape factor. 

The Penning source code uses the following conventions. For two electron 

excitation, the rate constant KPEN is set negative in calls involving the 

receptor specie. Bonor status is indicated by setting level index NEXP negative. 

In receptor calls, NEXP equals the total level Increment and is _<n* for one 

electron excitation/ionization, >n* for the two electron case in the receptor 

atom call. For the two electron case or single ionization, the receptor (single 

charge) ion type must also be called with the Penning effect invoked (LGP true) 

and with NEXP < n*. Inputs to the receptor ion call are obta ined from the atom 

call outputs. Consider the increment 6F^ in an appropriate element of vector 

F due to added Penning rates, ANMET contains for the atom at output. 

If u denotes the receptor ion terminal level (1 ^ u ^ » the receptor ,ion input 

ANMET = 6 6F . For ionization plus excitation, these are equal and the ANMET 
u u - 2 

value is simple passed in the calls; for single or double ionization, RPEN/N^ 

must be added to ANMET between calls. Also at receptor atom output (NEXP ^ n* 
at input), the' equilibrium density is multiplied by 3^^ and the ionization poten- 
tial is subtracted from the metastable energy. Both - these new forms are 
passed as inputs to the ion call. The final call to SRJM Involving the Penning 
effect is that of the donor. Inputs are: ANMET - the equilibrium donor ground 

level density based on the upper level density, ANMEQ - the actual ground level 
density, and RPEN - (B^fiFj^) for the donor which equals the negative of the same 
expression for the receptor (ANMET receptor atom output). During any call invok- 
ing the Penning effect, the excltation/ioniaation energy is subtracted from the 
EMET input. After the last receptor specie call, the increment in the electron 
rate heating is. thus the product of the output EMET and the rate RPEN. 
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-TFOR»SI SUB16 

SUBROUTINE SRJM { RL > JCB »F » AI »AR» ARR ♦QCRI »QCRR»KS»MS»AN»ANEP*QM»NPH, 
C ANHS*BETAPL.GLTE,CGME,D6ME»ZS,ZC,ANPL,LC,CVP»LGP,NC0RE,NEXP* 

C KPEN,RPEN»ANMET*ANMEQ»EMET»LERR» IGMX »FMD»CVOU-T ) 

INCLUDE SPEC1*LIST 

INCLUDE SPEC2»LIST 

INCLUDE SPECA*.LI.SJ 

INCLUDE SPEC6»LIST 

REAL JCB»KPEN,NPH(A) ,LPI 

LOGICAL LGP»LGI »LPRNT»LGL,LERR 

PARAMETER NWS=4«NLM»NPP=NPL+1 ♦KBV= (NPP*NPP) /4» 

C NSBT=2*NLM+1»NRF1=NPL+KBV+NLM 
DIMENSION RL(NLM*NPL) >JCB(NLM»NPL) *F (NLM) »ANX(NLM) ♦BETA(NLM) » 

C ANR(NLM) »LC(4) »BVS(KBV) »AM( NLMtNLM) *DW(NWS) »RMT(NLM»NPU » 

C CVP(2>4) ♦FMD(4) »CV0IJT(3»4) 

EQUI VALENCE (SV( 1) tANRCl) ) » (SV(NPL ) »BETA ( 1 ) ) » ( SV ( NSBT ) »BVS (1 ) ) » 

C (TMP(21)» DW( 1 ) ) ,(SV(NRF1) , AM{ 1* 1 ) ) , (TMP( 1 ) »RMT( 1»1) ) » 

C (SV(NSBT) ,ANX( 1) ) 

FX{X)=1.12*EXP{ -X)/(X+0,935*X**0.3) 

DATA S6/7,95775E-2/LPRNT/.TRUE. / 

Lll FORMAT( lH0,3X3iHIMPR0PER PENNING RECEPTOR LEVEL) 

L12 FORMAT (4X22HPENNING EFFECT* IGNORED) 

L13 FORMAT! 1H0,4X20HRATE MATRIX SINGULAR } 

114 FORMAT ( 1HO,4X24HJACOBIAN MATRIX SINGULAR) 

LERR=. false. 

IRST=MS 

LN=INX(KS»2)-1 

LG=INX{KS»5)-1 

LW=INX(KS»7)-1 

LF=INX{KS»8)-1 

L=INX(KS.l)+2 

GNS=DS(L) 

L=INDS(KS+1) 

LF=LF-L 

NB=NBAR(KS) 

• NL=NSTAR(KS)-1 

CALL CRATE(TE,KS,NL) 

DO 10 K=1»NL 
JCB(K»NPL)=RMT(K»NPL) 

DO 10 J=1,NL ' 

' RL(J»K)=0.0 

10 JCB( J»K)=RMT( J,K) 

DO 12 K=1»L 
IF(K.EQ.NB) MR=MS 
MS=MS+1 

12 ANR(K)=ANB{MS) 

ANPL=ANHS ' 

IF(ZC.LT.ZS) ANP,L=AN8(MS+1) 

S1=2.0708E-7/(GNS*TE»SQRT(TE) ) 

S2=15.79*ZC*ZC/TE 
S3=157.46*ZC**4*GNS 
S4=DS{NL+LN+1)»*2 
S4=S2/S4 ■ 

S5=1.0 

IF(ZC.LT.ZS) S5=ANEP*BETAPL 

AI=0,0 

AR=0.0 

ARR=0.0 72 

QCRI=0.0 
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QCRR=0.0 
DO 20 K=ltNL 
RSUM=0«0 
DO 22 J=1>NL 

RSUM=RSUM+JCB{K » J )*( ANR( K)-ANR( J) ) 

array F is used in the JACOBIAN MATRIX 
F ( K ) = ( S5*ANPL-ANR ( K ) ) *JCB { K»NPL ) -RSUM 
IFtK.LT.NB) F(K)=F(K)-ANR(K)*-DGME/ANEP 
Q=DS(K+LN)**2 
U=S2/Q 

A- S3»DS { K+LW ) *FX ( U-S4 ) / ( DS ( K+LG ) *0*DS ( K+LN ) ) 
Q=S1*DS(K+LG)*EXP(U) 

arr=arr+o*a 

ANX(K)=ANEP*Q*ANR{K) 

ARRAY BETA IS USED IN THE JACOBIAN MATRIX 
BETA(K)=Q 

A=A+ANEP*JCB(K»NPL) 

AR=AR+Q»A 

QCRR=QCRR+Q*U*A 

A=S5*A 

RL(K»NPL)=A 

SRMR=ANPL 

IF ( .NOT. { ZC.LT.^S) ) SRMR=SRMR+QM~ANC ( NCORE ) 

F(K)=F(K)+A*SRMR/ANEP 

A=Q*ANR(K)*JCB{K»NPU 

AI=AI+A 

QCRI=QCRI+U*(A+Q*RSUM) 

CONTINUE 

CALL LINWID(TE»T*ANEP,AN»KS,NL) 

CALL ESCAPE(ANX»KS,NL»AN) 

DO 25 K=l>4 
LC(K)=0 

IF( ICAVa»K).NE.KS) GO TO 25 
CVP( 1»K)=0,0 
CVP(2»K)=0*0 
DO 24 J=l»4 

IF( ICAV(2>K) .NE*LU(1»J) ) GO TO 24 
IF{ICAV(3»K).NE.LU(2,J) ) GO TO 24 
LC(K)=J 
CONTINUE 

IF(LC(K).EQ*0) LC(K)=-1 
CONTINUE 

S3=0,5101*S3/GNS 
DO 48 J=1»NL 

si=o.o 
$ 4 = 0.0 
RSUM=0.U 
U = 0.0 
SRMR=0.0 
R( J+IRST)=0.0 
R1=1.0/DS( J+LN)**2 
DO 46 K=1»NL 
IF(J.EQ.K) GO TO 46 
RL(J.K)=ANEP*JCB{ J,K) 

S1=S1+RL( J»K.) 

Q=R1-1.0/DS(K+LN)**2 

A=S3*Q*Q 

R2=DS(K+L6)/DS( J+LG) _73_ 

/rfta r Ji 
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R4=ANR(K)/ANR(J) 

IF(J.GT.K) GO TO 37 
A=A*ABS{DS(K+L*J+LF) ) 

IS=K+IFN( J*L) 
q=R3*A*PHI ( IS) 

RLI J»K)=RL< J*K)+Q 
L6I=(ABS(PHI(.I.S,)-1.*0-) .LT.l.OE-3 ) 
iF(LGI) AM( JtK)=Q*(PHI ( IS)-1,0) 

IFI.NOT.LGI) AM( J,K)=Q*W( IS) 

RSUM=RSUM+AM ( J , K ) *ANR ( K ) 

IF(.NOT.LGI) GO TO 27 
AMI J»K)=Q*WUS) 

SRMR=SRMR-R2*AM ( J *K ) /R3 

27 Q=1.0/(R3*R4)-1.0/R2 
IFCLGI ) 60 TO 28 

IF(ABS(Q) .LT.l.OE-7) Q=SI GN( 1*0E-7»Q ) 

SRMR=SRMR+AM ( j ♦ K ) / ( R3«^Q ) 

AM(J»K)=-AM( J,K)/{R2*Q) 

28 IF(NCAV.LE*0) GO TO 36 
DO 35 J1=1,NCAV 
I=LC(J1) 

IF(I.LE«0) GO TO 35 
IF(LU(1»I).NE.J) GO TO 35 
IF{LU(2»I ) .NE.K) GO TO 35 

CALL CAVITY{RC»RCP»GP(1,I ) »GP(2»I ) »GP ( 3» I ) »FMD( J1 ) ) 
CV0UT{1»J1)=GP(1»I) 

CV0UT{2»J1)=GP(2*I ) 

CVP(1*J1)=RC*NPH(J1) 

IF( (-R2*Q) •LE,l,0E-6) GO TO 29 
RCP=NPH{J1)*(RCP-RC) 

CVP{2»J1)=RCP 

RCP=RCP*R3/ANX(K) 

AM{ J»K)=AM( J,K)-RCP/{R2#Q) 

SRMR=SRMR+RCP/(R3*Q) 

RSUM=RSUM+CVP(2*J1 )/( ANEP*BETA( J) ) 
RC=R3*CVP(1,J1)/ANX(K) 

S1=S1-RC/(R3*Q) 

RL(J»K)=RL( J>tO-RC/(R2*Q) 

GO TO 30 

29 CVP{2»J1)=-1.0 

30 IGMX=K+IFN( 1,L) 

GP(3»I )=PHI (IGMX) 

35 CONTINUE 

36 A=S6*A 
IF(LGI) A=R2*A 
GO TO A5 

37 A=A*ABS(DS{ J+L*K+LF) ) 

IS=J+IFN(K,L) 

Q=A*PHI ( IS) 

S4=S4+Q 

L6I=(ABS(PHI( IS) -1.0) .LT.l.OE-3 ) 

IF(LGI) AM( J,K)=Q*(PHI(IS)-1.0) 

IFI.NOT.LGI) ami J»K)=Q*W( IS) 

U=U+AMU»K.) . 

IF{ .NOT.LGI ) 60 TO 38 
Q=Q*W( IS) 

AMI JtK)=R3*Q/R2 
SRMR=SRMR-Q _74_ 

38 Q=R3*R4/R2-1.0 
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ST=R2/R3-R4 
IF(LGI ) GO TO 39 

IF(ABS(Q) .LT.L.OE-7) Q=SIGN( 1.0E-7,Q ) 
IF(ABS{ST) .LT.l.OE-7) ST=SI6N( l.OE-7 ,ST ) 
SRMR=SRMR+AM( J*K)/Q 
AM( J*K)=AM( J»K) /ST 

39 IF(NCAV.LE.O) go TO 44 
DO 42 J1=1*NCAV 

I=LC( Jl) 

IFi I .EQ*0) 60 TO 42 
IFII.LT.O) GO TO 41 
IF(LU(1»I).NE.K) GO TO 42 
IF(LU(2»I).NE,J) go TO 42 
IFICVP(2»J1).LT.0.0) GO TO 40 
RCP=CVP(2*J1) 

RCP=RCP/ANX( J) 

AM( J*K)=AM( J,K)+RCP/ST 
SRMR=SRMR+RCP/Q 

RSUM=RSUM-CVP( 2 » Jl ) / ( ANEP*BETA( J ) ) 
RC=CVP{ l» J1)/ANX( J) 

Sl=Sl~RC/Q 

RL( J»K)=RL( J>K.)+RC/ST 

40 CVP(l»Jl)=CVP{l»JH-NPH(Jl)*FCAV 
CVP(2»J1)=A*ANX( J) 

GO TO 42 

41 IF( ICAV(2»J1).NE#K> GO TO 42 
IF( ICAV(3*J1).NE.J) GO TO 42 
GO TO 40 

42 CONTINUE 

44 A=S6*A 

45 IFtLGI ) A=A*wnS) 

IF(.N0T,L6I) A=A*(1.0”PHI(IS))/Q 

■ R( J+IRST)=R( J+IRST)+A 

46 CONTINUE 
S1=S1+ANEP»JCB{ J»NPL) 

RL( J*J)=-S1-S4 

AM( J»J)“SRMR 

F { J ) =F ( J ) + ( RSUM-U*ANR {'j 1 ) / ANEP 
R( J+IRST)=R( J+IRST)+S6*S1 
IF(J.GE.NB) GO TO 48 
Q=l,5+S2*Rl 

RL { J » J ) =RL ( J» J ) +0*GLTE-CGME’ 
Q=ANEP*(ZS+1*0-ZC)»ANR( J) 

DO 47 K=l»NL 

47 AM ( J » K ) =AM ( J »K ) -Q*BETA ( K ) * JCB( K.» NPL ) 

48 CONTINUE 

C FORM JACOBIAN MATRIX 

Sl=ANEP*ANEP»(ZS+1.0-Za/QM 
DO 49 J=1>NL 

JCB ( J » NPL ) =BETA ( J ) * JCB ( J » NPL ) 

DO 49 K=1»NL 

49 JCB(K»J) = RL(K»J )+AM(K»J>--Sl*F(K)*BETA( J) 

C START PENNING SOURCE CODE 

C PENNING EFFECT TEST 

IF(.NOT.LGP) GO TO 90 
LGI=«TRUE. 

IF(NEXP.GE.NSTAR(KS) ) GO TO 70 
IF(KPEN.LT,0,O.OR.NEXP.EQ.1) GO TO 80 
IF(NEXP.GE,NB) GO TO 60 
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50 IS=IABS(NEXP) 

R1=BETA{1)/BETA( IS) 

R2=ANMET*KPEN 
R3=ANMEQ*KPEN 
IF(NEXP.GT.O) GO TO 51 
R2=R2/R1 

R3=R3/R1 
NEXP--NEXP 
RA=-RPEN/BETA{ 1 ) 

60 TO 52 

51 RA=R2*ANR( 1)-R3*ANR(NEXP) 

RPEN=ANEP*BETA( 1) *RA 
RA=RA/ANEP 

52 RL(NEXP»1)=RL(NEXP»1)+R1*R2 
JCB(NEXP,1 )=JCB (NEXP,1 )+Rl*R2 
RL(1*NEXP)=RL( 1 »NEXP)+R3 
JCB( 1»NEXP)=JCB( 1»NEXP) +R3 

53 Q=1.0/DS( 1+LN) **2-1.0 /DS(NEXP+lN) **2 
F ( NEXP ) =F ( NEXP ) +R1*R4 

R{ NEXP+IRST)=R( NEXP+IRST)+S6*R1*R3 

RL(NEXP»NEXP)=RL{NEXP»NEXP)-R1*R3 

JCB(NEXP»NEXP)=JCBTNEXP,NEXP)-R1*R3 

54 LGL= (KPEN.GE.O.O.OR.NEXP.GT.NSTAR (KS ) ) .AND. (NEXP.NE. 1 J 
DO 56 J=1,NL 

S4=S1*R4*BETA( J) 

IF(LGI ) jCB(NEXP»J)=JCB(NEXPiJ) -R1*S4 
56 IF(LGL) JCB(1.J)=JCB(1,J)+S4 
IF{.N0T.LGL) GO TO 82 
RL{li>l)=RL( 1.D-R2 
JCB( 1*1)=JCB( 1,1)-R2 
F( 1)=F{ D-R4 

R( l+IRST)=R(l+IRST)+S6*R2 

EMET=EMET-S2*Q 

GO TO 90 

60 IFtLPRNT) WRITE(6»111) 

lprnt=. false. 

IFIKPEN.GE.O.O) GO TO 50 
LGP=. FALSE. 

WRlTE{6fll2) 

GO TO 90 

70 Q=1.0/DS( 1+LN) **2 

R2=ANMET*ABS(KPEN) 

R4=ABS(KPEN)/ANEP 

S3=2.0*ANPL 

IF{.N0T.(ZC.LT.2S) ) S3=S3-ANC(NG0RE ) 

S3=S5*S3 

IF(KPEN.LT.O.O) 60 TO 74 
72 S4= ANMET*ANR( 1 )-ANMEQ*S3 

RPEN=ANMET*ANR{ 1 )-ANMEQ*S5»ANPL 
GO TO 78 

74 IF(ZC.LT.ZS.0R.NEXP.EQ.NSTAR(KS ) ) GO TO 76 
KPEN=-KPEN 
NEXP=NSTAR(KS) 

GO TO 72 

76 IS=NEXP-NSTAR(KS) 

IFdS.LE.O) .60 TO 86 

IF { IS.EQ.NSTAR( K.S+1 ) ) GO TO 79 

IS=MS+IS 

S3=2.0*ANB{ IS) 
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RPEN=ANMET*ANR( 1)-S5*ANMEQ*ANB( IS) 

77 SA=ANMET«ANR( l)-S5*ANMEO*S3 

78 R4=S4*R4 
ANMET=RA*BETA(1) 

IF(NEXP.EQ.NSTAR(KS) ) R4=R4-ANMEQ*QM*A8S( KPEN ) /ANEP 
IF(NEXP.EQ. (NSTAR(KS)+NSTAR(KS+1) ) ) R4=R4-S5*ANMEQ*QM* 
C ABS(KPEN)/(2.0*ANEP) 

RPEN=ANEP*BETA( 1 )*RPEN»ABS{KPEN ) 

LGI=. false. 
anmeq=anmeo*bet A ( 1 ) 

GO TO 54 

79 IF(ZC.LT.ZS-l) GO TO 60 
S3=3.0*ANHS-ANC(NC0RE) 

RPEN=ANMET*ANR( 1)-S5*ANMEQ*ANHS 
GO TO 77 

80 Rlsl.O/BETAINEXP) 

81 R2=ANEP*BETA{ 1) 

R3=R2»ANMEQ*ABS ( KPEN ) 

R4=ANMET 

IFtNEXP.GE.NSTARIKS) ) GO TO 82 
GO TO 53 

82 EMET=EMET-S2*Q 
ANMET=R2 
ANMEQ=R1*ANMEQ 
GO TO 90 

86 R1=ANEP/S5 
GO TO 81 
90 CONTINUE 

C END PENNING SOURCE CODE 
J=l 

. IF(NCAV.lE.O) GO TO 96 
92 CONTINUE 

DO 94 K=L»NCAV 
I=LC(K) 

IF(I.NE.J) GO TO 94 
J=J+1 
GO TO 92 
94 CONTINUE 

96 IGMX=0 

97 IF(J.6E.,5) GO TO 100 
IF(GP(1»J) .GT.O.O) GO TO 98 
GO TO 100 

98 IF( I6MX.EQ.0) IGMX=J 
IS=LU( l.J) 

IS=LU(2»J)+IFN( IS»L) 

GP(3»J)=PHI (IS) 

J=J+l 
GO TO 97 

100 IF(NB.GE.NSTAR(KS) ) GO TO 160 
IS=NB-1 

NV=NSTAR(KS)-NB 

CALL INVERS(RL»8VS,NV».TRUE.tF,$130) 

DO 120 K=1.NV 

ANB(K+MR)=-RL(K+IS*NPL)*ANPL 
DO 120 J=1,IS 

ANB(K+MR)=ANB(K+MR)-RL(K+IS»J)*ANB{ J+IRST) 

F( J)=F( J)-RL( J»K+IS)*F(K+IS) 

RL( J»NPL)=RL( J»NPL)-RL( J»K+IS)*RL(K+IS»NPL) 

DO 120 I=1*IS 
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120 RL( J*I )»RL( J*I )-RL( J,K+IS)*RL(K+IS*I ) 

60 TO 135 
13CJ WRITE(6»113) 

133 LERR=«TRUE, 

GO TO 160 ’ 

135 CALL INvERS(JCBfBVS,NV, .FALSE. .ANR»$150; 

DO 140 J=1,IS 
ANR{ J)=BETA{ J) 

DO 140 K=lfNV 

ANR { J ) = ANR ( J ) -BETA ( K+ I S ) * JCB ( K+ I S ♦ J ) 

* JCB( J*NPL)=JCB( J»NPL)-JCB(K+IS.NPl)*JCB(»C+IS» J) 

DO 140 1=1, IS 

140 JCB( I*J)=JCB( I,J)-JCB(I ,K+IS)*JCB(IC+IS,J) 

GO TO 160 
150 WRITE(6*114) 

GO TO 133 

160 return 

C IF NBAR IS LESS THAN NSTAR THE F0Li_0WIN6 SUBROUTINE' IS USED IN 

C CALCULATING EFFECTS OF QUASI-ST ATIONARY LEVELS 

SUBROUT I NE I NVERS ( AMS » B VS » NV, LOGC » F 
LOGICAL LOGC 

REAL AMS{NLM,1) »BVS(NV,1) ,F(NLM) 

DO 10 K=1»NV 
DO 10 J=1,NV 

10 AMI J»K)=AMS( J+IS.K+IS) 

DO 20 J=1,NV 
DO 14 K=1,IS 

14 BVSI J»K)=AMS( J+IS»K) 

IFILOGC) BVS(J,N8)=AMS( J+IS»NPL) 

20 IFILOGC) BVSI J»NB+1)=FIJ+IS) 

I = IS 

IFILOGC) 1=1+2 

CALL SOR I AM *NLM » N V ♦ B VS ♦ NV » I ♦ $30 » DIn ) 

GO TO 40 
3 ^ RETURN 6 
40 DO 50 J=1,NV 
DO 44 K=1*IS 

44 AMSI J+IS»K)=BVSIJ»K) 

IFILOGC) AMSI J+IS,NPl)=SVSI J,NB) 

50 IFILOGC) FI J+IS)=BVSI J,NB+1) 

RETURN 

C THE FOLLOWING SUBROUTINE CALCULATES FREQUENCY AVERAGED CAVITY 
C RATE FUNCTIONS NEGLECTING HOLE BURNING 
SUBROUTINE CAVI T Y I RC ,RCP »G»X,P»CF ) 

CF=CGM*P 

CK=I 1.0-LPI IX) ) /I l.O+CF) 

S=-G/ALOGlCMR) 

ST=AMAX1I 1.0-S»0.0) 

RC=FCAV*S*I 1.0-CK*ST**0.3 
ST=AMIN1(S,0.99) 

IFIS.LE.0.99) go to 20 
IFIS.GE.1.01) GO TO 10 
CK=CK*EXPt 1.0E3*IST-S) ) 

GO TO 20 
1^ CIC=0.0 

20 RCP = FCAV*S*I1.0-CK*I1.0-1.j>» 5I ) /il.u-i.! M 

RETURN 
END 
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Subroutine MATRIX is used to collect the results from calls to SRJM and form 
the overall rates and Jacobian* The argument list is as follows: 


ANR 

vector of dependent variables (lower level densities, photon 
densities, dilnX^/d^, input 

FF 

vector of lower level rates 

DMJ 

Jacobian matrix 

BETA 

vector of lower level Saha factors 

ANEP 

predicted input 

ANP 

vector of maximum ionization densities 

GLTE,CGME,DGME 

(see SRJM list). 

DE 

electron rate heating/kT^ 

LCD 

input integer vector; total number of lower level equations 
of species with indices <i is LCD(i) 

AN,EP,DI,DEL 

(see CONSRV list) 

LERR,CVP 

(see SRJM list) 

ANPH 

photon density vector 

FMD 

vector of mode spacing parameters 

CVT 

same as CVOUT in SRJM (cavity outputs) 

GOUT, LOUT 

maximum positive noncavity gain outputs: GOUT(l,k) is gain 
logarithm, G0UT(2,k) is escape factor, L0UT(i,k) are level 


indices (i = 1 lower, i = 2 upper) for ion or atonj types k 


In the subroutine DEM is the overall rate matrix* Also RPEN is the Penning 
rate constant rather than KPEN and PRATE is the net Penning rate rather than RPEN 
(see SRJM list) • DPEN is the increment to the electron rate heating due to the 
Penning effect. At output, ANEP is the corrected electron density. 
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SUBROUT INE MATR IX ( ANR »FF » DMJ»BE TA , A.NEP ♦ ANP » GLTE »CGME »DG^'E OE» LCD ♦ 

C AN»EP»DI»DEL*LERR»CVP»ANPH»FMD»CVT,C0L'T,LCUT) 

INCLUDE SPECltLIST 
include SPEC2*LIST 
INCLUDE SPEC4,LI5T 
INCLUDE SPEC5»LIST 
include SFEC6»LIST 
INCLUDE SPEC?* LIST 

PARAMETER NLL=NCM+1 ♦NSTl = 3*KM^NPL.NST2=NSTl + KM*KM»f\S7 3 = NST2 + 2*NH30 

PARAMETER N17S=XM+ 1 ,.N17S2 = N173+KM ,N1 ?S3 = N1752 + KM 

dimension ANP(NCM) »LCD('iLL)»S^Tr(NLM} »D£M(KM»KM) tD.v.j (KM»KK} » 

C FF (K.M) »2ETA(KM) »GGCK.K} *HlKy) >BTP I KMI'tRL (.\LM»NPLJ » 

C AJ5(.\LM»\PL) »F ,LC(4) »C\/P{2»4) tFMD(4) *CVT13»4) » 

C GOUT ( 2 » iNSP ) » LCUT I 2 » NGP ) 

LOGICAL LERR.LGP 


lQJ I VALENCE C BE TP ( 1 ) , SV ( 1 ) ) » ( GG { 1 ) »3TRMTX ll ) ) » (H( 1} »STRMTx(iM7S) ) 
C (3TF(1)»£TRMTX(N1702) ) » ( Ftll tSTRKTX (N17S3 ) ) » (DEM(1»1)» 

C STRMTX CMSTl ) ) » ( RL( 1 » 1 ) »£TRMTX( NST2 } ) » ( AJB( 1»1 ) »STRMTX ( NST3) ) 

MS=C 


> 


is=o 


DC 5 J=1»NC 
iriJ.EQ.NC) GO TO 1 
KL=IC{ J+l)-l 
GC TC 2 


1 X2=NSP 

2 K1=IC{J) 

• CO 5 K.^K1»k2 
L“INuw(K+l) 

MI=NBAR IK)*1 
DO 3 I=1»MI 

3 ANC( I+Mi£)^ANR(l4 IS) 

I S- 1 S+M. I 

KS=MS+L 
5 CONTINUE 

DO 7 K=1»IS 
DC 7 J=1»IS 
D£IUJ.K)=C.O 
7 DMJ(J»K)-L.O 

CALL CCNSRV(DPEN.ANP,EP,DI ,DEL.GM,AN ,C) 

B 1 = 0 » G 

B2-0.C 

B3=L.0 

DPEN=C,0 

R1=0,0 

R2=0.0 

LR1 = C 

LR2 = G 

LD1 = 0 

LD2-C 

MS=C 

IS =0 

DO ICO J=:,NC 
IFiJ.EQ.NC; GO TO 31 
K2-MU J + D-l 
GO TO 51 
31 K2=NSP 
51 ia=IC{J) 

25^K2-R1+1 
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ANHS=ANP( J) 

NEXP=C) 

L1=LCD( J)+l 
L2=LCD( J+1) 

DO ICO K=K1,K2 
ZC=K-K1+1 

IF(K.LT.K2) 60 TO 61 
BETAPL=1.C/ANEP 
GO TO 71 
61 L=INDS(K+1) 

‘ BETAPL=BT{MS+L+1) 

71 LGP=. FALSE, 

L=NBAR(K)-1 
DO 81 I=1,L 
BETA(I+IS)=BTII+MS) 

81 GG{ I+I5)=ANEP«ANEP*(ZS-2C+1.0)/QM 
IF(.NOT.LPN) GO TO ^0 
PENNING SOURCE CODE ~ SET A 
IPEN=C 

IF( J.EQ.JRC.AND.K.EQ.Kl) GO TO 10 
IF{ J.EQ,JRC.AND*K,EQ,K1+1) GO TO 20 
IF( J,EQ.JDC.AND.K,EQ,K1) GO TO 30 
60 TO 40 
10 IPEN=l 

RATIN=RPEN 

LGP=.TRUE, 

nexp=jrl 

i F ( RPEN,LT. 0. 0 ) NEXP=NEXP+NSTAR ( K) 

L=IC(JOC) 

I=INX(L,2)-1 

EMCT=15.79/TE 

EMET=EMET/0S( I + 1)»»2-EMET/D3t I + JDl.)»»2 

L=l.r 1 . 

MI=0 

DO 12 I=1»L 
12 MI=MI+INDS( I+l) 

ANMET=ANEP»BT ( M 1+ JDL ) *ANB { iMI+JDL ) 

ANMEQ= ANEP^f BT ( M 1 + JDL ) *ANB ( MI + 1 ) 

A1=ANEP*BT(MS+1) 

ANGR=Ai*ANB{MS+l) 

IF(NEXP,LT,NSTAR(K) ) ANUR=Al»tANB (MS+ JRU 
IF(NEXP.EQ.NSTAR(K) ,AND,K,EQ.K2 ) ANUR=A1»ANHS 
GO TO 4C 
20 IPEN=2 

RATIN=RPEN . 

lF{NEXP*Lf,NSTAR(Kl) ) GO TO 4 
LGP=,TRUE, 

NEXP=1 

IF(RPEN.LT,O.OJ 60 TO 22 
ANUR=A1*ANEP*BT ( MS+1 ) *ANPL 

GO J0_ 40 ... _ 

22 NEXP=JRL 

IF(NEXP,EQ,NSTAR(K) ) 60 TO 24 

ANUR=A1*ANB(MS+JRU 

GO TO 40 

24 ANUR=A1*ANHS 

GO_TO 40 

30 IPEN=3 

RATIN=ABS(RPEN) 
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PRATE=SAVE 

NEXP=-JDl 

EMET-w* G 

IF (RPEiN.lt. 0.0) a.N<JR=ANMET*ANUR 

ANM£T = AiNSjR 

ANMEa=ANGR 

C END SOURCE SET A 

40 CONTINUE 

CALL SRJM(RL»AJS»F»AI »AR» ARR.3CRI »QCRR.!C»MS»AN»ANEP,C:M,ANrH. 

C ANHS'.BETAPL .GlTZ,CQME»DSME»ZS,2C ♦ ANPL ♦ LC »CVF »LoP » J *'iNEXP .RAT I N . 

C PRATE » ANMET » ANN'EO . EMET » LERR » IGMX ♦ FMD > CVT ) 

IF(LERR) GO TO IIC 

GOUTa»lC)=G.O 

GOUT(2>K)=0.0 

LOUT(i»K)=0 

LCUT(2»K.)=0 

IFUGMX.LE.C) go to 41 
G0UT(1.K:)=GP(1» IGMX) 

G0UT(2»K)=GP( 3. I6MX) 

L0UT{1»K)=LU(1.IGV.X) 

L0UT12»K)=LU(2. IGVX) 

41 A2=K2+1-K 

A3 = ANEP»BETAPL*A.\PL 
B1=B1+A2*(ANEP«AI-A3*AR) . 

32=B2+A2^^A3*( ARR-AR) 
rF(K.EO.K2) 52=32+AR«ANC( J) 

33 = B3+A3*GCRR-ANEF.^QCRI 
L=NBAR(K)-1 
00 42 MI=1,L 
BTP(MI+IS)=BETP (MI ) 

H{MI+IS)=-F(MI ) 

IF(K.EQ.K2) H(MI+IS)=H(MI+I3)+QM«RL(MI »NPL) /a 

42 GG(MI+IS)=GG(MI+IS)*BETP(MI ) 

DO 50 MI=1*l 

IF(<.LT.K2) GO TO 46 
FF(MI + IS)=RL(MI .NPLX-ANP( J) 
do 44 NI=Ll»L2 

. 44 DMJ(MI+IS»NI)=-RL(MI,NPL)*ANEP*BTP(NI ) 

GO TO 48 

46 FF(MI+IS)=O.C 

DEM(MI+IS»L+1+I3)=RL(MI .NFL) 

DMJ(MI+IS»L+1+IS) =RL(MI »NPL) 

48 DO 50 NI=1.L 

DEM(NI+IS»MI+I5i=RL(NI .Mil 

DMJ(NI+IS»MI + I3) =DMJ( NI + IS.MI + IS)+AJS(N1 ».MI > 

50 CONTINUE 

NI=LCD(NC+1) 

A3*A2*ANEP 
A2=A2*AR 
DO 60 1=1. NI 
IF(K.LT.K2) GO TO 54 
DO 52 MI=L1»L2 

52 DMJ( I »MI ) =DMJ( I .MI )-ANR( I ) *A2*AN£P*B TP { MI ) 

• 60 TO 56 

54 DMJ( I.L+1+I3)=DMJ( I.L+1+IS)+A2«ANEP*dETAPL*ANR(I ) 

56 IF( I.GT.IS.AND. I .LE.IS+L) GO TO 6C 
DO 58 Ml=l.L 

DMJ( I»MI-H3)=CMU( I.MI + IS)-ANR(I ) »A3* A JB<MI »NPL) 
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6C CONTINUE 

' IF(,N0T*LPN) go to 99 
C" PENNING SOURCE CODE - SET 3 

IF(IPEN.GT.O) GO TO t 70 » 80 »90 ) » I PEN 
GO TO 98 

70 IF(LGP) GO TO 72 
LPN=*FALSE. 

GO TO 99 
■ 72 a2=0.0 

A3=0.0 
AA=0 »0 

SAVE=PRAT£/ANEP**2 

LR1=IS+1 

LR2=IS+JRL 

IF(RPEN.LT.O.O) LR2=lR2+L 

IF(RPEN.GE,0*0*AND.JRL.GE,NBAR{ K1 ) *0R .RPEN. LT. 0* O.AND. 

C JRL.GE*NBAR(K1+1 > ) LR2=0 
IF(NEXP.GE.NSTAR(Kln 60 TO 73 
R1=SAVE 
ANMET=SAVE 
60 TO 76 

73 R1=ANMET 

IF(RPEN.LT*«J.O.AND.NEXP.LT.NSTAR(Kl).+NSTAR(!U+ln GO TO 74 

IF{RPEN.GE.G.0.AND.!<;i.*EQ#K2) A2=RPEN 

IF(RPEN.LE.O.O) A2=ANEP*BETAPL*ABS(RPEN)/2.0 

A2=A2*ANMEQ*QM/ANEP 

A3=SAVE 

74 SAVE=ANMET-A2 
ANMET=ANMET+A3 

76 R2=ANMET 

IF( JRL.EO.NSTAR (Kl + 1) . AND.RPEN. LT . 0. 0 .OR* JRL .EQ. NSTAR ( Kl ) .AND. 

C K1.EQ.K2) 60 TO 78 

GO TO 79 . 

78 R2=0.0 

A4=ANMEQ*ANEP»ABS(RPEN) 

IF {RPEN.lt. 0.0) A4=A4*ANEP*SETAPL 

79 IF(RPEN.LT.0.0.0R.K1.LT.K2) GO TO 93 

80 DPEN=EMET*PRATE 
60 TO 98 

90 LD1=IS+1 

LD2=IS+JDL 

IF{ JDL.GE.NBARIKl) ) lD2=0 
GO TO 98 

96 DO 97 I=L1»L2 

IF( J.EQ.JRC) DMJ{ LR1» I ) =DMJ(LR1 ♦ I )+A4*DTP( I )/BETA(LRl) 

IF{ J.EQ. JDC ) DMJ{ LDl, 1 ) =DMJ(LD1 ♦ I )-A4*BTP (1 ) /BETACLDl ) 

IF( J.EQ*JDC.AND.LD2.GT.0)DMJ{LD2»I)=DMJ(LD2,I )+A4*BTP( I )/3cTA(LD2) 

97 CONTINUE 
GO TO 99 

98 IF(K.EQ.K2.AND.NEXP.N£.0) GO TO 96 
C END SOURCE SET B 

99 IS=IS+L 
100 CONTINUE 

IS=0 

- NI=LCD(NC+1) 

DO 1C4 K=ltNSP 
L=NBAR{K)-1 
DO 102 MI=1»L 
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A3=0.0 

IF(J*EQ.lR 1) A3=R1/BETA{LR1) 
IF{J.EQ.LR2) A3=-R2/BETAaR2) 
IF(J.EQ.LDI) A3=-R1/BETA{LD1) 
IF{J.EQ.LD2) A3=R1/BETA(LD2) 

DO 102 1=1, NI 

IF{I.GT*IS.AND. I.LE.IS+L) GO TO 102 
DMJ( J»I ) = DMJ( J, I ) + (H.( J-)-+A3 )*G6( I ) 

102- CONTINUE 
104 IS=IS+L 

CALL C0NSRV(ANEP,ANP,EP,D1 »DEL»QM,AN, 1 ) 

CGME=B1 

DGME=81+B2 

DE=B3*ANEP+0PEN 

DO 106 J=1»NI 

DO 106 K=1,NI 

106 FF( J)=FF( J)+DEM( J,K)*ANR(K) 

110 RETURN 
END 
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Function TMSA is used to calculate the spontaneous emission factor for the 
cavity. (DM is the mode spacing parameter.) 

Subroutine EICP calculates electron-ion collision parameters according to 
the results of Ref. 26, (see also Refs. 1 and 5.) The ratio of collision 
integrals B* Is not explicitly used in the program and is. set equal to a nominal 
value of unity. Argument Z is the ion charge number. 

Subroutine ELEN calculates dAnT /dt at a given time TAU (see argument lists 

e 

of SRJM, MATRIX), The electron energy equation is discussed in Refs. 1 and 5. 

Input forcing functions are given by FORIN in the following manner; output BM 

is the magnetic field intensity (nominally null), IDIS is an integer that defines 

the type of forcing function, FT is T and FP is d£nT /dt if IDIS ^0, FT is the 

total electric fielc tangent to the magnetic field and FP is the perpendicular 

electric field if IDIS » 1 and FT, FP are the component current densities if 

IDIS ^ 1, The dummy subroutine calls ADUMYn are to be replaced by subroutines 

gi_ving electron— atom collision parameters (elastic cross section Q and collision 

Integral ratios B* and C*). Nominal (final) variables in TEMPS common are: 

H, HI, H3 - Hall parameters (see references), SIGMA - electrical conductivity, 

ALPHT - thermal diffusion ratio, FEC - electron-heavy parti '.le elastic energy 

8 ^3 “1 

exchange collision frequency, Q - joule dissipation (units of 10 Jm sec ) 

and QR — continuum radiation (MKS units if multiplied by 5.4 ^ 10 ), 
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-TFOR»SI SUB18 

FUNCTION TMSA(ETA»PS,DM,CL»CAR.WL) 

TMSA=ETA 

RETURN 

END 

-TFOR_»^ SUB19 ■ 

SUBROUTINE EICP ( TE»ANE,Z»QEI »BEI »CEI ) 
BEI=1.C 

DL= -864».88*SQR-T-( T-E/ANE ) 
TS=0.151*DL*TE/Z 
. IF(TS-.LE.l.O) GO TO 1C 
DL=DL/TS 

QEI=DL*DL*ALOG( 3.0»TS) 

TS-TS**0,2 

CEI=0.351*TS/{TS-0,325) 

RETURN 

10 QEI=1.0986«DL«0L/TS 

CEI=0.52 
RETURN 
END 
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SUBROUTINE ELEN ( GLTE ♦ TAU» DE. COME »AN» ANP ) 

INCLUDE SPEC1,LIST 
include SPEC2»LI5T 

DIMENSION CSTR(8) *A ( 2 . 2 ) »F ( 2 ) »T N ( 2 ) » AN ( NSP ) ,ANP{NCM) 

EQUIVALENCE (A( 1 ♦ 1 ) »CSTR ( 1 ) ) , ( T N ( 1 ) ».CSTR I 5 5 ) ,(F{1) jCSTRtT) ) 
EQUIVALENCE ( TMP( 21) »CSTR( 1) ) » ( TMP(29 > »S) , ( TMPC 30) .Q) » ( IMP ( 31) »B) 
C {TMP(32)»0 »(TMP(23)»ALPHT) » ( TMP { 34 ) ♦ DEL ) » ( TMP ( 35 ) ♦HI), 

C (TMP(36) ,H3 ) , (TMPC 37) ,H) , (TMP(38 ) ,SI6MA) , ( TM? ( 39 ) ♦ FEC ) , 

C (TMP{40),QR) ,(TMP(41) ,Z) ♦ (TMP ( 42 ) ♦ AM) 

CALL FORINITAU, IDIS,BM,FT,FP) 

IFdOIS.LE.O) GO TO 20 
H=15*79/TE 
GR=0.0 
FEC=0.0 
DO 1 J=l,8 
1 CSTR(J)=O.C 
DO 10 J=1,NC 
IF(J.EQ.NC) GO TO 2 
K2=IC( J+l)~l 


GO TO 3 

2 K2=NSP 

3 Kl = IO(J) 

Hl=K2-ia 
DO 8 K=K1,K2 
ZsK-Kl 

I=INX(K»1)+1 
AM=DS(I ) 

C=DS( 1+2) 

IF(K*LT.K2) S=AN(K+1) 
IF(K.EQ.K2) S=ANP(J) 

I = INX(K,2)+NSTAR(i<)-l 
Q=(Z.+1.0)/(DS( I )-0.5) 

■ QR=QR+S*C*CZ+1.0)**4*EXP(H*Q*Q 
S=AN(K) 


IF(K*GT.K1) GO. TO 6 
1 = 1 

C DUMMY SUBR CALLS ARE TO BE CONVERTED BY COLLECTOR 
JFLJ-.EQ.D CALL ADUMYl ( TE»Q»B»C ) 

IF(j.EQ.2) CALL ADUMY2 ( TE,Q,B,C) 

IFU.EG.3) CALL ADUMY3 ( TE ,Q,B,C ) 

IFIJ.EQ.4) CALL ADUMY4 ( TE »Q,B»C ) 

1FIJ.EQ.5) CALL ADUMY5 ( TE ,Q,B,C ) 

IF(J.EQ.6) CALL ADUMY6 ( TE ,Q,B,C ) 

aO-JTCLj . 


6 ' 1 = 2 

CALL EICP(TE,ANE,Z,Q,B,C) 
. 7 IN.m=TN(I)+S 
S=Q*S 


F-{ n =F { I ) +S 



A( 1»I)=A(1,I)+B*S , 


A(2»I )=A(2,I )+C*S 


. JLELIA8^(Hlr2).*CT,0.5} i 


2*H1+1*0 


S=ANRU) 


GO TQ 6 

. 8 

CONTINUE 

, 10 

CONTINUE 


GO TO 8 
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DO 12 K=l>2 
S=F{K) 

DO 12 J=l»2 
12 AC J»K)=A( J,K)/S 
S=FC2)/F( 1) 

Q=1.0+S 

B=1*0+0.7*S-0.48»A(1,1) 

C=i.2*CA(2,l)+S*A{2>2) )-Q . 

ALPHT=C/B 

DEL=ALPHT*C/Q 

H1=SQRT(TE) 

GR=ANE*QR/H1 
FEC=FEC*H1 
H1=H1*(F( 1)+F(2 ) ) 

SIGMA=3.402E3-»ANE/ril 
Hl=2. 1233E3»BM/H1 
H3=0.4«Hl»Q/8 
Q=DEL/( 1,0+H3*H3) 

H=H1+Q*H3 

Q=1,0-Q 

H=H/Q 

DEL=1.0-DEL 

IF( IDIS.EG.l)' G=SIGMA*( FT5^FT/(Q*( 1,0+H»H) )+FP*FP/OEL) 

IF( IDIS.GT.l) Q=(Q*FT*FT+DEL*FP»FP) /SIGMA 

GLTE=DE/( i.5*ANE)+(4.829E-3*Q-2.6085E-5#GR)/(ANE«TE)-9,088E-4» 
C FEC*(TE-T)/TE-CGME 
GO TO 30 
'20 TE=FT 

GLTE=FP 
30 return 
END 
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